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ABSTRACT 


This report describes models and computer codes which may be used to 
describe flow reactors in which high purity, solar grade silicon is produced 
via reduction of gaseous silicon halides. A prominent example of the type of 
process which may be studied using the codes developed in this program is the 
SiCli*/Na reactor currently being developed by the Westinghv^use Electric Corp. 

During this program two large computer codes were developed. The first 
is the CHEMPART code, an axisymmetric, marching code which treats two-phase 
flows with models describing detailed gas-phase chemical kinetics, particle 
formation, and particle growth. This code, based on the AeroChem LAPP (Low 
Altitude Plume Program) code can be used to describe flow reactors in which 
reactants mix, react, and form a particulate phase. Detailed radial gas-phase 
composition, temperature, velocity, and particle size distribution profiles 
are computed. Also, deposition of heat, momentum, and mass (either particulate 
or vapor) on reactor walls is described. The second code is a modified version 
of the GENMIX boundary layer code which is used to compute rates of heat, 
momentum, and mass transfer to the reactor walls. This code lacks the 
detailed chemical kinetics and particle handling features of the CHEMPART code 
but has the virtue of running much more rapidly than CHEMPART, while treating 
the phenomena occurring in the boundary layer in more detail than can be 
afforded using CHEMPART. 

These two codes have been used in this program to predict particle forma- 
tion characteristics and wall collection efficiencies for SiCl/, /Na flow 
reactors. It is found that large input enthalpies (large H-atom inputs) are 
required to prevent Si(^) droplet formation. (This enthalpy is supplied by 
introducing large quantities of arc-heated hydrogen in the Westinghouse 
reactor.) On the other hand, large hydrogen flows mean short transit times of 
gas through the reactor and hence short times for wall collection of Si. It 
is anticipated that an important application of these codes will be their use 
in finding operating conditions where droplet formation may be minimized and 
high collection efficiencies may still be realized in reactors of the 
Westinghouse type. 


NOMENCLATURE! SECTIONS II-V 


damping constant in Van Driest 's mixing length formula 
wall skin friction coefficient 

mass diffusivity of species i (Fickian for vapor, Brownian for 
particles) 

number of monomer units (approximate measure of particle size) ; 

see Eq. (24), Section II. 

mass flux of species i (kg m“^ s) 

^thermal diffusion ratios” defined in Eqs. (32) and (33), 

Section II. 

pseudo "specific rate constant” defined by Eq. (45), Section IV 
(due to thermophoresis) 

Knudsen number of species i (= T/R^) 

mean free path of gas; also used for turbulence mixing length 
molecular weight of "carrier gas" 

molecular weight of species i 

molecular weight of diffusing vapor (i.e., silicon monomer) 

Nusjelt number 
pressure 

radial distance measured from reactor tube centerline 
radius of ith species (molecular or particle radius) 

Reynolds number 

Schmidt number (= y/5 D^) of species i 

Stanton number 

temperature 

streamwise component of gas velocity in the developing boundary 
layer (do\mstream reactor region) 

radial component of gas velocity in the developing boundary layer 

(doimstrerm reactor region) 

thermal settling velocity (Section II) 

streamwise distance along reactor tube length 
mole fraction of species i 
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Greek 




‘^TTl 




e . /k 
1 

X.,A 

e 


M 

V 

P 


Pi 


a. 

X 


do\<mstream distance measured from reactor tube v/all 
mass fraction of species i 


fraction of molecules undergoing "diffuse reflection" from particle 
surface. 

thermal diffusion of thermophereti'' factor for species i 

momentum accommodation coefficient 

value of in the continuum limit (Kn^ = 0) 

value of in the f ree-molecule limit (Kn^ = 

exchange coefficient for heat transport 

effective exchange coefficient for mass transport 

species and thermal boundary layer thicknesses 

eddy diffusivity 

energy-well-depth parameter in Lennard-Jones potential law 

thermal conductivity of particle and carrier gas 

dynamic viscosity of gas 
kinematic viscosity of gas 
density of gas 
density of particle 

diameter of particle of species i 


Subscripts 

b bulk property 

c centerline of carrier gas 

CE Chapman- Enskog 

eff effective 

h pertaining to heat transfer 

i species 

m pertaining to mass transfer; also used to indicate conditions at 

edge of species layer 
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mp 

0 
r 

SEM 

W 

1 


melting point 

initial condition 

radial or reference quantity 

Stokes-Einstein-Millikan 

quantity evaluated at the reactor tube wall 
referring to monomer 


TP-392 


ai/2 

A 

b 1 / 2 
B 





AG 
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NOMENCLATURE! SECTIONS VI-VIII * 

defined by Eq, (71), Section VI. 

constant in expression for (see Eq. (61)) (RC(K,1)) 

defined under Eq. (63), Section VI. 

activation energy (see Eq. (61)) (RC(K,3)) 

specific heat of mixture (2 X c^ )(CPBAR(K)) 

i i i 

specific heat of ith species (CPTB(I,K)) 

diffusion (Brownian) coefficient for particles, ith mass class 
turbulence energy dissipation (TURDIS) 

ratio of actual drag coefficient to Stokes flow drag coefficient 
for particles, ith mass class 

moles of ith species per gram of fluid (ALPHA(I,K), RALPHA(I,K)) 
number of particles per gram of fluid, ith mass class 

Gibbs xree energy of ith species at standard state (1 atm) 
(GTB(I,K)) 

ratio of actual heat transfer coefficient to that for Stokes flow, 
ith mass class particles 

change in standard Gibbs free energy for a reaction, 2 (gj_) 

products - 2 (gj^) reactants ^ 

i 

Gibbs free energy of formation 
particle thermal velocity (GPP(I,J)) 

static enthalpy of mixture (H, STATEN) 
enthalpy of ith species (HTB(I,K)) 

heat of formation of ith species at T = 298 K (HE (I)) 
stagnation enthalpy of mixture (HSTAG, STAGEN) 


* Where a FORTRAN variable in the CHEMPART code exists corresponding exactly 
to the tabulated quantity, the FORTRAN name is given in parentheses. 
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K 



Kn. 
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Le 

m, 
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Mi/2 

N 

NPG 


Nu 


P 

Pr 

r 

r . 
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r 

in 

ri/2 

TC 1 f 2 

R 

R 

w 

Re 

Re 

s 


St 

T 



enthalpy of formation 
Boltzmann’s constant 

fon^ard rate coefficient (RATE) 

backward rate coefficient 

eddy viscosity coefficient 

eddy viscosity coefficient for Donalds on/ Gray Model (see Eq. (70)) 
coagulation rate coefficient (RATEP(I,J)) 

Knudsen number of particles, ith mass class (XKNUN(I)) 

equilibrium constant (K) 

latent heat of vaporization 

Lewis number (laminar or turbulent) (XLE(K)) 

mass of particles, ith mass class (MASSP(I)) 

Mach number at half radius, defined under Eq. (70) (QQ300) 
temperature exponent in reaction rate equation (Eq. (61)) (RC(K,2)) 
number of particle mass classes 
Avogadro’s number (AV) 

particle Nusselt number 
static pressure 

Prandtl number (laminar or turbulent) (SIGMA(K)) 

radial coordinate normal to jet or reactor centerline (Y(K)) 

particle radius, ith mass class (RZ(I)) 

inner mixing zone radius (QQ200) 

defined under Eq. (66b) (QQIOO) 
defined under Eq. (65) (QQIOO) 
universal gas constant (R) 
reactor radius (RWALL) 

Reynolds number 

particle Reynolds number, ith mass class (REP) 

supersaturation ratio 
particle Si.okes number 
static temperature (T(I),RT(I)) 

particle temperature for ith mass class (UP(I,K)) 
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Greek 

a 



Y 

6 . 

X 


e 

£ 


n 

X 

y 

y. 


V 

p 


wall temperature (TWALL) 

X component of velocity (U(K),RU(K)) 

X component of particle velocity, ith mass class (UP(I,K)) 
r component of velocity 

rate of production, of ith species (WDOT(I)) 

rate of particle production, ith mass class (QXP(I,K)) 

rate of production from jth reaction 

molecular weight of mixture (2 F.)”^ (inverse of WTMIX(K) 

i 

molecular weight of ith species (WTMOLE(I)) 

coordinate parallel to jet centerline (X) 
mole fraction of ith species 

distance from wall, (r — r) 

mass fraction of ith species 


constant for external control of eddy viscosity (see Eqs, (63) to 
(69) (XK2) 

particle capture efficiency 
surface tension 

particle mean free path, ith mass class 

eddy diffusivity for turbulent flow; defined as y/p 
emissivity (EPS) 
defined by Eq. (63) 
molecular mean free path 

effective viscosity for turbulent flow (XMU(K)) 
molecular viscosity (UM(K)) 

eddy viscosity 

kinematic viscosity 
gas density (RHOG(K)) 
particle density (RHS) 

Stef an-Boltzmann constant 


ix 


TP-392 




stream function (PSI) 


Subscripts 


e 

i 

j 


Pi 

o 

w 


evaluated at edge of mixing layer (free stream) 
ith species 

value at nozzle (jet) exit 
particle mass class i 
evaluated at axis of symmetry, r = 0 
evaluated at wall 


Miscellaneous 


( ’O, 


2 

i 

[ ] 


absolute value 

partial derivative with respect to P; y being held constant 
summation over i ii’pecies 
concentration of gaseous species 
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I , INTRODUCTION 

This report covers the development of models and associated computer codes 
x^hich describe processes and phenomena involved in novel methods for the pro- 
duction of high purity silicon, A prominent example^ of such a process is one 
undergoing development by the Westinghouse Electric Corp, in which the produc- 
tion of silicon is performed via a chemical reaction between gaseous silicon 
tetrachloride and sodium jet streams, carried out in an arc-heated, turbulent 
flox^ tube reactor. In this process the reaction betx\reen SiCl/* and Na is 
completed in an upstream reactor section and silicon vapor and/or droplet 
formation occurs. The process then strives for efficient separation and 
collection of the silicon in a downstream section by the collection of silicon 
on the reactor x^alls. In order to provide a reasonably complete description 
of processes such as this the following computer codes were developed or adapted: 

GENMIX^ - Solves the developing, turbulent boundary layer problem 
encountered in the downstream reactor section due to the flew of silicon- 
containing hot product gases through this cooled-wall portion. Deposition of 
silicon on the walls is taken to be controlled by convective-diffusive process- 
es a Thus separation/collection of silicon vapor can be described by this code. 
The code accounts for variable fluid properties but neglects radiation and 
thermophoretic (Soret) effects in computing the reactor’s vapor collection 
efficiency. 

MPDEU - Solves the particle transport equation in the downstream reactor 
section. This equation describes the mass transfer rate of silicon droplets 
controlled by convection, Fick and eddy diffusion, as well as the Soret effects. 
Since the gas-phase velocity and temperature fields in the developing tubular 
flox^ of the reactor govern the abovementioned transport processes, GENMIX is 
coupled to MPDEU. 

CHEMPART - Solves the turbulent, jet mixing problem including detailed 
gas-phase chemistry and particle formation. This code also treats mass, momen- 
tum, and energy flux to the walls. This code is an extension of the well-knoxm 
LAPP^ code for rocket plume exhaust studies. 

The development of these codes is based on a rationale best explained x^^ith 
the help of Figs. 1 and 2. Figure 1 shoxv^s an SiCl^/Na flow reactor with a 
simple reactant input geometry. The significant phenomena occurring x^ithin 
the reactor are specified, along with rough values for the temperatures prevail- 
ing in each region. Specifically, the following phenomena must be described if 
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the fate of the silicon produced in the reactor is to be computed; (i) in the 
upstream mixing region, rate of mixing and reaction, (ii) also in the upstream 
region, rate of Si droplet formation, (iii) in the do^mstreara, developing flow 
region, rate of boundary layer growth, (iv) also in the do^mstream region, 
rate of particle growth via agglomeration and heterogeneous condensation, and 
(v) also in the domstream region, rate of diffusion and thermophoretic trans- 
port of Si vapor and droplets to the reactor walls o 

In order to describe the situation of Fig, 1, the scheme indicated in 
Fig, 2 has been employed. In this schematic, three regions are indicated i In 
the first, upstream region in which the jet of Na vapor is mixed with the silicon 
halide, the CHEMPART code is used to extract detailed information concerning 
the mixing and rate of reaction. Particle formation, if supersaturation occurs, 
is also computed. In the next section, where boundary layer development occurs, 
CHEMPART is again used to follow particle growth rates, slow gas-phase chemistry, 
and wall deposition. However, due to the size and expense of running this code, 
it was felt that another smaller code would be desirable. Hence, a boundary 
layer code x^as developed, based on the GENMIX code, x^hich could treat the 
boundary layer problems without detailed chemistry, but x^ith a mass transport 
model to compute x^zall deposition rates. Finally, some 50 or more diameters 
doxmstream, a fully-developed flox^7 occurs. At this point, one-dimensional 
calculations could be used to compute Xv^all deposition rates. CHEMPART tests 
for this situation and once a uniform core flow develops, detailed calculations 
are performed only at the centerline and xs^ithin the boundary layer. Although 
not in fact a one-dimensional calculation, this technique results in much the 
same result in terms of economy. (GENMIX runs rapidly enough that such modifi- 
cation is unnecessary.) 

In the sections x^hich follox^, the theoretical background needed for the 
use of these codes is described. Then the code themselves are explained. 

Finally, their use in studying SiCl^/Na flow reactors is described. Listings 
of the codes and examples of input and output are given in Appendices. 
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II . MODELING OF THE REACTOR BOU ND A R FLAYER 

A major problem faced in the production of silicon via gas-phase reaction 
of silicon halide in flow reactors is the collection of the Si vapor or drop- 
lets formed by the reaction* One means of collecting the silicon is to allow 
enough residence time within the reactor for diffusion of the silicon to 
reactor walls followed by condensation of vapor (or capture of droplets) on 
the walls. Such a condensation method is being employed by Westinghouse. ^ 

In view of the potential significance of this condensation process, a 
major effort in this program has been the adaptation of a powerful boundary 
layer code to permit the calculation of rates of transport of Si through the 
growing boundary layers of such flow reactors. In this section the theory of 
mass transport through such boundary layers is described. Since it is envi- 
sioned that flow reactors that depend on Si vapor wall deposition^ will operate 
above the Si dew point, emphasis will be on the calculation of vapor transport. 
However, since droplet formation may occur in the boundary layer or in the 
doxmstream part of the core flow within the reactor, the transport of a 
particulate phase is also fully examined. 

A, THE BOUNDARY LAYER MODEL 

The basic tool used for study of mass transport through the reactor 
wall boundary layer is the well-knoxsm GENMIX code.^ This code is used to 
obtain the information concerning the velocity and temperature distributions 
near the reactor wall needed to compute rates of mass transport through the 
boundary layer. The code is applied to the do^^stream section of the reactor 
where it is assumed that chemical reaction is complete and the gases are thus 
effectively inert. The physical problem thus involves an analysis of a hot, 
turbulent developing flow field within a circular pipe with cooled walls. Due 
to the large temperature variation between the reactor core flow and the walls, 
it is essential to include the effects of non-uniform fluid properties and the 
possible formation of silicon droplets due to condensation (caused by cooling 
silicon vapor). Hox^ever, for simplicity, at present we restrict attention to 
the case with no condensation. That is, silicon is assumed to exist only in 
the vapor state. 
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Yet another simplification is rendered in this analysis by viewing the 
overall product gas mixture as an effective binary system, with silicon as 
one species and all other gases (mainly NaCl, Ar, H 2 ) as the only other species 
(referred to subsequently as the "carrier" gas). The molecular weight and 
other properties of the carrier gas are taken to be a weighted average of the 
properties of the individual components. On this basis, the diffusional be- 
havior of silicon vapor through the surrounding gas can be uniquely described 
in terms of the appropriate mass dif fusivities (or mass transport coefficients) . 
On the other hand., the thermal diffusivity (or heat transport coefficient) of 
the overall product gas mixture will be nearly equal to that of the carrier gas 
alone, since silicon concentrations are typically low. In the case of silicon 
vapor deposition, the relative importance of heat and mass transfer processes 
inside a typical reactor is nearly equal. That is, the Lewis number for the 
system is close to unity. The effect of turbulence, as treated in the model 
considered, lies mainly in enhancing the diffusive transport of heat, mass, 
and momentum.* The enhancement due to the so-called turbulent counterparts 
of the abovementioned dif fusivities will be assumed to be such that the 
"effective Lex^is number" remains close to unity. That is, turbulence enhances 
heat and mass transfer rates equally. However, according to current belief 
and experimental data the enhancement of momentum transfer rates due to tur- 
bulence is somewhat less (a fact represented by setting the turbulent Prandtl 
number equal to 0.9). 

It is clear, therefore, that a good description of turbulent transport 
coefficients is pivotal to the success of the present model. For now, a mix- 
ing length approach to turbulence modeling has been adopted. The chief drax^- 
back of this approach is that it relies heavily on an empirical determination 
of several constants, some of xijhich might vary from one flow situation to an- 
other. This handicap will have to* be overcome by the use of more sophisticated 
(differential) models of turbulence. This latter class of turbulence models, 
x^^hich describe the mixing length via differential equations (rather than alge- 


The diffusive transport of momentum (or vorticity) is analogously controlled 
by the kinematic viscosity (v), or momentum diffusivity. 

The Prandtl number is the ratio of momentum diffusivity to thermal 
diffusivity. 
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braic ones, as done here), can, of course, be blended into the present formula- 
tion. However, we believe that since turbulent heat/mass transfer predictions 
using even the most elaborate approaches sometimes incur errors of up to 15%,^ 
retaining the mixing length model would at least give us the advantage of great- 
er simplicity with relatively little loss of achievable accuracy. Moreover, 
the mixing length model has been extensively studied and tested in pipe flow 
configurations against both experimental data^ and the predictions of more 
sophisticated models,® 

In view of the abovementioned assumptions and recognizing that the develop- 
ing pipe flow processes of interest are characterized by comparable changes in 
the radial and axial directions (i.e., tv;o-dimensional, axis3rTmnetric) one may 
adopt the Patankar and Spalding^ turbulent boundary layer formulation for the 
conservation of overall mass, momentum, energy (or stagnation enthalpy) and 
silicon vapor species.* No attempt will be made to reproduce these well- 
established parabolic, partial differential equations here since they have 
been fully detailed in Refs. 2 and 5. Instead here we choose to discuss 
aspects not immediately obvious from their discussion. 

Implicit in the governing equations used by these authors^ are the 
following assumptions : 

(i) The densitv fluctuations caused by turbulence are of insignificant 
importance compared to the flux contributions of other fluctuating quantities 
(e.g., velocity, temperature, concentration). 

(ii) The dissipation rate of turbulent kinetic energy (due to velocity 
fluctuations interacting with the mean velocity) is negligible compared to 
the viscous dissipation rate of mean kinetic energy. 

Under these assumptions, it is possible to reduce the otherwise complicated 
turbulent boundary layer equations to a form similar to the better-understood 
steady, compressible, laminar boundary layer equations, provided one defines 
suitable "effective” values for the momentum, thermal, and mass dif fusivities 


It can be estimated from the magnitudes of the various dif fusivities that 
the velocity, temperature, and species boundary layers will all develop at 
nearly the same rate. The "fully-developed" state is not reached until 
some 40-60 diam doxmstream from the pipe inlet. 
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discussed earlier. These have been defined by the following flux-gradient re- 
lationships (assuming Newtonian behavior)'^: 


8u 

^ ~ ^eff 3y 

a" = -r C 

^ h,eff p 8y 

9Y. 

i " = -r — - 

■^i i,eff 9y 


(shear stress) 

(1) 

(heat flux) 

(2) 


(mass flux of species i) (3) 


where is the effective viscosity, the ratio of effective thermal 

conductivity to the constant pressure specific heat of the carrier gas 

(C^) , and the product of the mean gas density (p) and the effective 

mass diffusivity (D^ of silicon vapor in the carrier gas. u, T, and 

denote the time-averaged streamwise velocity, absolute temperature, and mass 
fraction of silicon, respectively. 

Each of the effective transport coefficients above is the sum of a 
laminar (or molecular) contribution and a pseudo, turbulence- induced con- 
tribution. The latter may be greater by some two orders of magnitude in 
fully turbulent (i.e., nearly inviscid) regions, such as the reactor core 
flow. However, closer to the reactor walls turbulence energy is rapidly de- 
pleted due to the dissipating influence of molecular viscosity, and transport 
is affected mainly by the laminar (molecular) mechanism. As a result the 
turbulent parts of the transport coefficients can vary significantly across 
the reactor cross-section. The motivation for introducing a "mixing-length,” 
or any other turbulence model, is to describe this variation realistically. 

In the present analysis the following radial variation of mixing length, 
from the reactor axis (or centerline) to the wall, was prescribed: 


A6 


Ai 

K 


I = 


Ky[l-exp(-y"''/Ab ] , 0 ^ y ^ ^ 


(core layer region) 
(near wall region) (4) 


* We neglect for the moment (see the next section) thermophoretic transport. 
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where y is the distance from the wall, 6 is the thickness of the turbulent 

-ft 

boundary layer, and X, K, and A are taken to have the constant values 0.09, 
0.435, and 26.0, respectively. Physically, such a distribution of mixing 
length distinguishes regions of high local Reynolds number, y = yu^/v,‘ (or 
nearly inviscid flow) from those regions near the wall which are governed by 
the progressively increasing effect of molecular viscosity (i.e., a "damping” 
of the mixing length with decreasing y"^, according to Van Driest 's’'® exponen- 
tial law) . Note that this permits a viscous sublayer region next to the reac- 
tor walls so that the influence of even weak turbulent fluctuations is felt. 
Just outside the viscous sublayer the mixing length obeys an undamped, "defect 
law" behavior (i.e., fi, = Ky) . 

Using the above mixing length distribution it is possible to express the 
effective viscosity as: 


eff 


^turbulent 


laminar 


(5) 




9u [ 

3yl 


+ y 


( 6 ) 


where |— | is the magnitude of the time-averaged velocity gradient normal to 
dy 

the reactor wall. Since this quantity decreases away from the wall (becoming 
zero at the boundary layer edge) turbulence is operative only in regions in 
which the velocity gradient exceeds a small critical value. If, in any flow 
region the turbulence velocity scale | £ falls below a certain fraction of 
the local velocity, this quantity is set to the local prevailing velocity. 
Such a circumstance tends to occur in the fully-turbulent core flow region 
of the reactor. 


It has been shown®’® more recently that A. is a sensitive function of pressure 
gradient and wall blowing or suction. This parameter determines the thick- 
ness of the "viscous sublayer" region. 

" u^ is the so-called "friction velocity" defined here in terms of the local 
shear stress by the relation u^.^ 5 (r/p)^/^ 
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B. COMMENTS ON THE NUMERICAL SOLUTION PROCEDURE 

Since the governing equations to be solved are parabolic partial differ- 
ential equations, their solution requires the specification of appropriate 
initial and boundary conditions. In the present study the initial conditions 
at the pipe entrance were specified as uniform (corresponding to centerline 
values) profiles of u, T, and across the pipe cross-section. The domain 
of integration was taken to be the region between the pipe centerline and 
the wall. Boundary conditions on u, T, and were specified at the wall as: 

U = 0 (7) 

T = T (constant) = 1700 K (8) 

w 

Y. = == 0* (9) 

X 

The last condition is a consequence of assuming that silicon vapor is in 
equilibrium at the wall (i.e., the partial pressure equals the saturated 
vapor pressure) . The centerline boundary conditions on these variables were 
taken to correspond to the symmetry conditions: 


9u 

ay 


0 , 


ay 


0 , 


aY 

ay 


0 


( 10 ) 


Using the above initial and boundary conditions, GENMIX produced solutions 
at specific downstream stations along the reactor by marching forward in steps 
(whose size was proportional to the local boundary layer thickness) . At each 
step the solution to the coupled system of conservation equations is carried 
out using a fully implicit, 6-point, finite difference algorithm (based on an 
integral approach that ensures conservation, rather than the usual Taylor se- 
ries expansion) . The algorithm uses a substitution method for solving the re- 
sulting system of coupled algebraic equations (with a tri-diagonal matrix) , 
thus avoiding the time-consuming and unreliable operation of inverting the 
coefficient matrix. It was found in this study, as in others, that GENMIX is 


A and B are constants with the values 46710.0 K and 7.3166 x 10^° N m 
respectively. 
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extremely efficient as a code. For instance, the full solution for a 9 m 
reactor length with some 80 to 100 radial gridpoints is completed in about 
17 s (CP time on CDC 7600 computer). 

The remarkable efficiency of GENMIX is really a result of the many dif- 
ferent time-saving features built into the code. While the program chooses 
its own forward marching step size, the cross-stream grid spacing is specified 
as input. In the present problem, the gridpoints were so chosen as to be 
unequally spaced over the pipe radius, with a greater resolution capability 
near the wall, where gradients were steep. Furthermore, it must be noted that 
the solution in GESMIX is not carried out in physical space (x,y) but rather 
in the transformed von Mises^^ coordinates (x,m) . In this latter coordinate 
system, the normalized stream function m always provides a fixed integration 
domain: 0 ^ cu 1. In general, the advantage of working in this transformed 

space is the greater solution accuracy and ease that results from imposing 
boundary conditions at fixed extremities. However, in the solution of con- 
stant radii pipe flows (such as the present case), this advantage is not fully 
realized since the two extremities remain fixed even in physical space. Of 
course, an alternative solution tactic may have been to solve the problem be- 
tween the wall and the edge of the developing boundary layer. In this latter 
case, working in the transformed stream function ordinates may be a definite 
advantage. 

A major disadvantage of the GENMIX solution procedure may be its handling 
of the near-wall flow. The reactor wall is a point of singularity in the 
(x,oj) space because the flow velocitv there is zero and it multiplies the 
highest derivatives of the transformed equations. This disconcerting feature 
also tends to render the transformation process back to physical space somewhat 
inaccurate gridpoints close to the wall. In order to circumvent some of 

the difficulties associated with this singular point, the Patankar-Spalding pro- 
cedure generates solutions for the region next to the wall using a turbulent 
Couette flow analysis (further simplified using a constant property assumption) . 
Thus one does not obtain a rigorous finite difference solution of the two- 
dimensional boundary layer equations in the region close to the wall. However, 
this device greatly enhances the numerical solution efficiency since wider grid 
spacings can be used. Moreover, by arranging the region of Couette flow very 
close to the wall, the error due to this approximation may be rendered acceptable. 
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C. SILICON TRANSPORT LAWS 

Using GENMIX a description of the boundary layer flow field can be 
generated. The transport of Si to the walls nay then be generated once the 
diffusion constant describing its Fick and Soret transport properties are 
obtained. 

The mechanism of transport of silicon particles (including, as the 
smallest particles, vapor molecules) onto reactor walls is dependent on the 
size of the particles. For particles larger than about a micron, deposition 
is controlled by turbulence- induced inertial impaction; sub-micron sized 
silicon droplets deposit by a convective-diffusion process. In order to pre- 
dict the mass flux to the reactor wall of the latter class of droplets it is 
necessary to recognize two separate contributions to the time-averaged diffu- 
sional flux: (i) a concentration gradient-driven (Fick) flux and (ii) a 

temperature gradient-driven (Soret) flux. The relative importance of these 
contributions itself varies with particle size, the former being dominant for the 
smaller sub-micron sizes while the latter controls the behavior of particles with 
sizes closer to a micron, as sho\^ qualitatively in Fig. 3. Thus, in general, 12-1 a 


78-142 



FIGURE 3 VARIATION OF DEPOSITION 'RATE OF SILICON PARTICLES WITH SIZE 
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for a species i (vapor or particle) that is dilute in a carrier gas the dif- 
fusional mass flux (per unit area, per unit time) at any point, jl^, (correspond- 
ing to a particular pressure and temperature) is given 





^ [(VY^) + (1 - Y^) I (T£) 


pD 


+ j 


turb 


( 11 ) 


where p is the local mass density of the carrier gas, the local mass fraction 
of the diffusion species, T the local absolute temperature (with the operator V 
denoting gradient of the particular scalar, temperature or concentration) . The 
first two terms in Eq. (11) represent the dif fusional mass fluxes due to concen- 
tration and temperature gradients, respectively. The third term, j'turb’ 
mass flux contribution due to turbulent or ”eddy diffusion." Assuming, as is 
customary, that the turbulent diffusion also has a Fickian nature, one can write, 


' turb 


-PCp(7Vp 


( 12 ) 


where is the so-called "eddy dif fusivity . " Particles larger than about a 
micron would inevitably be deposited on the walls due to inertial impaction 
caused by turbulent velocity fluctuations normal to the wall. Hence ^ for such 

“ 4 '* 

large particles the overall mass flux to the wall would be dominated by j^urb* 
However, the sub-micron size particles will be driven close to the wall by tur- 
bulence eddies before the concentration and temperature gradient driven fluxes 
(i.e., Pick and Soret effects) become important in the final stages of deposi- 
tion, within a nearly laminar (viscous) sublayer next to the wall. 

However, before deposition rate (or mass flux) can be calculated, one 
needs to establish the Pick (or Brownian) diffusion coefficient D^, the dimen- 
sionless thermal (Soret) diffusion or "thermophoretic" factor a^, and the "eddy" 
dif fusivity The former two transport pr’operty coefficients depend upon the 

molecular nature of the species present in a given "effective" binary mixture 
(e.g., silicon-argon, in the case considered later in this report), and on the 
prevailing pressure and temperature. In addition, they are both strong func- 


Under typical silicon reactor operating conditions the concentration 
dependence of these transport coefficients will be negligible, due to the 
diluteness of silicon in the carrier gas. 


12 


r«*i 


TP-3: 2 


tions of particle size, changing by several orders of niagnitude (as shown later) 
in traversing the particle size spectrur/i from molecular diameters to about a 
micron. 

Interestingly enough, while the values of the Fick dif fusivities (i.e., 

D. and e ) are always positive, ensuring that Bro^mian and turbulent transport 

ip 

of Si always occur in the direction of decreasing concentration (towards the 
reactor vzalls) , the value of ot^ may sometimes be negative. The latter circum- 
stance is usually characteristic of diffusing species that are lighter than 
the carrier ga^ (e.g., silicon vapor (28.09) diffusing through argon (39.94)). 

Thus the Soret flux has the potential to drive silicon vapor towards higher 
temperatures, away from the "cool" walls, and thereby reduce the overall col- 
lection efficiency of Si. Fortunately, as shown later, the value of a^. for 
Si particles (defined here as any molecular cluster exceeding roughly twice 
the mass of a single vapor molecule) is invariably positive, since the parti- 
cles are heavier than the surrounding gas. That is, while condensed Si drop- 
lets will always be driven towards the reactor walls, uncondensed Si vapor 
could actually be driven away from the walls due to the Soret effect. 

Therefore, in addition to considering the mechanisms of Broraian and 
turbulent transport of silicon, we investigate the conditions under which Soret, 
or temperature gradient-driven diffusion, can be important. Although thxs 
latter effect has been ignored in the past in the belief that its contribution 
to the overall mass flux would be negligible, we conclude that the Soret effect 
can indeed play a crucial role in the Si (particle/vapor) separation process, 
under typical reactor conditions. Fortuitously, the Soret flux for condensed 
Si droplets will always favor deposition onto the walls. 

For the sake of comparison, it is necessary to define the "eddy" diffusivity, 

£ , for mass transport. In general, by viewing turbulence as the random agitated 
p 

movement of small packets of fluid, it is possible to def ine ip terms of a 

"mixing length" (in much the same spirit as the "viscous mean-free-path" in the 
kinetic theory of gases represents an average distance for effective transfer 
of momentum by molecules). Thus, 


Particle size is described in terms of diameter or radius since it will 
be assumed here that the particles (silicon liquid droplets) are perfectly 
spherical. 
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Cp ^ u'il. -V 1^1 (13) 

where )li is a mixing length, u' the large scale turbulent velocity fluctua- 
tions present in the core flow of the reactor, and j-^l the magnitude of the 
time-averaged (steady) normal velocity gradient • Both h and 1^| will vary 
with distance above a wall, y, with li being zero at the wall since turbulence 
is precluded (i.e., the ”no-slip” condition). I'Tliile the many available models 
of turbulence provide different expressions for the variation of h with dis- 
tance normal to a wall, accurate estimates of and the proportionality constant 
in the above expression, for any given flow situation, usually result from ex- 
periments. Lin et al^^suggest the following expression, valid in the viscous 
sublayer flow, next to a wall: 


, < 5 


where v is the kinematic viscosity and y is a non-dimensional normal distance 
from the wall defined as: 


y'*’ = |y U(f/2)'^=]/v (15) 

with U the average velocity and f the Fanning friction factor. For particles 
small enough to follow the small-scale eddy motions close to a wall, the eddy 
diffusivity does not depend upon particle size. 

Although silicon particles of widely varying sizes must be considered in 
the modeling of the deposition processes within any proposed reactor, reliable 
information about the "molecular" transport coefficients (i.e., and a^) is 
available only at the two extremes of the size spectrum. For particle sizes 
much smaller than the mean-free-path of the carrier gas, the well-kno\m 
kinetic theory results of Chapman and Enskog^^ (CE) apply. In this limit of 
"f ree-molecule" flow, surrounding gas molecules do not suffer appi’eciable 
changes in their distribution functions upon collision with the particle. 
Furthermore, since expressions for are an outcome of CE theory using 
second-order terms of the perturbation series, while the expressions come 
out of first-order terms, one concludes that will be more sensitive to 
interactions between colliding molecules. Thus more sophisticated molecular 
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potential models than the Lennard-Jones 12:6 (e.g., the exponential model)^^’^® 
might ultimately be necessary to accurately determine in this limit, al- 
though the Lennard-Jones might suffice for engineering purposes. 

For particles much larger than the mean-free-path of the surrounding gas, 
the gas structure becomes "invisible” and the behavior can be described by the 
Navier-Stokes equations for a "continuum,” subject to the usual "no-slip" 
boundary conditions at the particle surface (e.g., the Stokes-Einstein (SE) 
expression for particle diffusivity^^) . 

The main difficulty in specifying transport properties is encountered 
for particles in the intermediate size range, on the order of a mean-free- 
path. This presents the notorious ".usolved problem of the "transition regime." 
Although several attempts have been made to extend the continuum results into 
the transition regime, by imposing "slip" or jump- type boundary conditions at 
the particle surface, these are at best heuristic. Understanding of the real 
transition mechanisms involved is altogether insufficient. However, the transi- 
tion regime experimental data of Millikan^ ° and others have made it possible to 
establish various empirical correction factors (as discussed later) that cover 
a wide portion of the transition regime. Unfortunately, even such widely used 
interpolation formulae as the Stokes-Einstein-Millikan (SEM) expression are 
unsatisfactory since they do not exactly match the predictions of CE theory 
when extended to the "f ree-raolecule" limit. 

Crucial to the accurate modeling of silicon deposition, therefore, is 
a proper description of the abovementioned transport coefficients. Since 
present knowledge regarding the "transition regime" is insufficient, it would 
be justifiable to suggest universal interpolation formulae for these coefficients, 
covering the entire particle size (or Knudsen number*) range. Moreover, 
mathematically continuous functions that blend together the known results in 
different regions would be both physically and computationally more desirable. 


The particle Knudsen number (Kn^^) will be defined here as the ratio of the 
mean-free-path (£) in the surrounding gas to the particle radius (R^) (i.e., 
Z 

Kui = -^) . Thus Knj^ 0 corresponds to the "continuum" limit while Knj[ 

refers to the "free-molecule" limit. It is sometimes assumed^ ^ that the 
range 0 < Knj_ < 0.25 represents the region of "slip-flow" while 0.25 < Kni < 
10 is the region of "transition-flow. " 
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In what follows, a novel rational approacl*^^ to creating such "universal” 
descriptions of the Pick diffusivity (i.e., D^) and the Soret factor (i.e,, 
a^) is outlined. 

1. Universal Formula for 

The expression for diffusivity adopted in this report blends the near- 
continuum SEN relation (valid for small Knj[ values) with the free-molecule CE 
expression (valid for large values of Kn^, > 10) according to the following 
universally valid relationship^ ® : 


D. 

1 


D 


SEM 


+ 


CE 


[1 -f exp {- 


~ ^SEM 

C(Kn^ - Kno)}] 


(16) 


x^here C and Kno are constants, 


the transition from to as the Kn^« is increased. 

oET'l CjE 


The magnitude of C controls the abruptness of 

Kno is the Knudsen 


number value at the "point of inflection" of D 


SEM' 


Figure 4 illustrates the 


rationale behind choosing such an interpolation function to describe the dif- 
fusivity of particles in a gas. 
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16 



TP-392 


The Stokes-Einstein-Millikan expression for dif fusivity , D__„, is based 
2 0 

on Millikan^ s famous oil-drop experiments which permitted measurement of 
the isothermal drag force on particles of various sizes moving through a gas. 
Using a stochastic analysis of Bro^mian motion, Einstein^^ had earlier related 
the particle diffusivity to the "mobility,” or relative velocity per unit force, 

through the reaction: 


D. = kTB (17)* 

1 


the mobility, B, being given by the Stokes formula for lox^7 Reynolds number 
continuum flow around a sphere: 


B 


1 


(18) 


where y is the dynamic viscosity of the gas. Recognizing that as the particle 
size approaches the gas mean-f ree-path, the drag for a given velocity becomes 
less (due to "slip") than predicted by Stokes law (i.e., the mobility increases), 
Millikan suggested the following expression for mobility, which is a well- 
accepted correction of the Stokes result: 


B 


1 

6'rryRi3 


(18a) 


where 


B = [1 + C Kn. + C^Kn. exp(- 0.88/Kn.)]"" 

m X d X X 


(18b) 



(18c) 

(18d) 


with a denoting the so-called momentum accommodation coefficient, a measure 
m 

of the efficiency of momentum transfer between the particle and the surrounding 


5V 

k is the Boltzmann constant. 
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gas molecules. Combining Eqs. C17) through (18d) yields the needed expression 
for diffusivity, valid at least for small departures from Kn^ = 0: 


D 


SEM 


kT 


[1 + C Kn. + C.Kn. exp(- 0.88/Kn.)] 
m 1 d X 1 


(19) 


At the other extreme, for particles much smaller than the gas mean-free- 
path (but not necessarily as small as molecules) the diffusivity should be 
given by the Chapman-Enskog relation^ 


D 


CE 


1.8583 X 10““^ 




( 20 ) 


where p is the pressure (atm) , and the molecular weights of the carrier 
gas, c, and diffusing molecular species, i, respectively, the effective 

collision diameter (angstroms) and Q the reduced collision integral. Equation 
(20) is strictly valid only for molecules with spherically symmetric potential 
force fields, with the supplementary relations; 


a 


ci 


1 

2 


(Oc + a.) 


a 



(See Ref. 23 for a tabulation 
of this function.) 


( 21 ) 


e . 

Cl 






where a . and e , can be taken to represent the size and energy-well depth in. 

Cl Cl 

say, the Lennard-Jones molecular interaction model. 

To generalize Eq. (20) for the case of small particles, we need to redefine 
only the molecular weight and the molecular diameter Treating the par- 

ticle as a "heavy molecule," g times the mass of the corresponding molecule, 
its "molecular weight" can be expressed as: 


Values of am depend on the nature of the particle and the surrounding gas, 

19 

Experimentally determined values for various substances usually lie in the 
range 0.9 < cim < 1.0. In this report a^ = 1 (perfect accommodation) is 
assumed for the results presented later. 
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= g Ml = 



( 22 )* 


where Mi is the molecular weight of the substance constituting the particle, 

Ml 


g is the ratio of the particle volume to the molecular volume (i.e., 


of the substance, and is the particle density, 
diameter of a particle is simply 


The "molecular” 


Vi 


) 


a. = (2R.) X 10 (A) 

1 1 


(23) 


if the particle radius is expressed in centimeters. Combining Eqs. (20) 
through (23) yields the following expression for the diffusivity of very 
small particles: 


D 


CE 


1.8583 X 10 



(24) 


where 


Cl 


1 

4 


+ 2 X 10® R^)' 


(24a) 


0. ^ 



(24b) 


g 



(24c) 


Incorporating the expressions for and given by Eqs. (19) and (24), 

into Eq. (16) yields the required, universally valid, equation for the diffu- 
sivity D^. The value of Kno was taken to be 0.44 so that the point of inflec- 
tion in Eq . (16) coincided x^ith the inflection point of the SEM formula (i.e., 
Eq.(19)). The value of C was chosen to be 1.5 for the results presented below. 


is the Avogadro number. 
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2. Universal Formula for 

Using a rationale similar to the one adopted for but necessarily 
different because entirely different (and more complicated) expressions 
usually describe we propose the following universally valid formula 
for the Soret factor: 


“i = ( l A-Kn T)°» 

X X 

where aoo is the "thermal diffusion factor" predicted by CE theory for Kn^ 
and ao is the value of the Soret factor as Kn^ ->• 0. The constant A is deter- 
mined by forcing the prediction of Eq. (25), in the transition regime, to 
agree with established experimental data for some value of Kn^ of order unity. 

It is known that a particle suspended in a gas with a uniform tempera- 
ture gradient will drift through the gas under the action of a thermal force. 
\^en this "thermophoretic" force just balances the opposing drag force on the 
particle, a constant thermal settling velocity is achieved. Using Millikan ^s 
oil-drop apparatus, modified to provide for a temperature gradient, Phillips 
was able to accurately describe the thermal settling velocity up to Kn^ « 0(1), 
for particles of widely varying thermal conductivities (including the coupling 
that usually exists between the velocity and temperature fields surrounding a 
moving particle). We have incorporated^^ the Phillips formulae into a computer 

• 4 - 

program that predicts the thermal settling velocity (V^) some value of the 
Knudsen number (say, Kn^ = 1) and used this to find the corresponding Soret 
factor (cip). Knowing aco and ao, as shown below, A can then be determined as: 

A = — ° - (26) 

a=o - “P 

for a given particle-carrier gas combination and given flow conditions. 

Since the CE expression for is extremely complicated, we utilize the 

analysis of Waldmann and Schmitt which yields the following expression for 
->• 

Vp valid for small particles whose drift in a temperature gradient can be 
associated with the stronger bombardment by higher energy molecules from the 
"hot side": 
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a .D . 


X X 


^ '’(l+ScJ 


X 


(27) 


where is the diffusivity already discussed, v the kinematic viscosity, and 

a, the fraction of molecules that undergo ”dif fused reflection” from the par- 
d 

tide surface. It turns out that the above result for does not exactly 

match the CE result for molecular sized particles. In order to achieve this 

match a new constant correction factor a is introduced as follows: 

corr 


a 



(1 + 1 


1 


a 

corr 


(28) 


Physically, the effect of a is to account for the variation in a, with 
increasing particle size.’' The ratio (v/D^) is the particle Schmidt number 
(Sc.). It is a measure of the relative importance of momentum to mass trans- 
fer. Since the diffusivity of large particles is small, unlike gas molecules, 
Sc^ is typically large (as seen later. Table I). 

Here ao is taken to be related to the thermophoretic velocity deduced 
by Epstein, for particles much larger than the gas mean-free-path (using 
a slip-flow procedure). Using Epstein^s result one can write: 


Uo 



1 + 


1 

0.5(A./A ) 
X c 


(29) 


where and A^ are the thermal conductivities of the particle and gas, 
respectively. 

Equations (26), (28), and (29) when incorporated into Eq. (25) provide 
the required universal formula for the Soret factor a^. 


3 . Transport Properties 

The analysis just presented for the evaluation of the transport 
properties of particles dilute in a carrier gas, was first applied to the 
case of silicon particles in argon, under typical reactor operating condi- 


V’Then (a.)^„ is negative, a ^1 may be assumed. 
X C>iii corr 
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tions. Figures 5 and 6 show the numerically computed values of and 
over a wide range of Knudsen number (^j^) • Table I is provided to gain an 
idea of the corresponding particle size in terms of the other possible size 
variables: g, or Sc^. Note that while g and are absolute measures 

of particle size, Kn^ and Sc^ depend strongly upon the prevailing pressure 
and temperature. Table II gives the values of input constants and the com- 
puted values of some useful quantities. The following features regarding 
Figs. 5 and 6 are also noteworthy: 

(i) Both and change rapidly by orders of magnitude, in going 
from very small (g = 1, = 1.58 A, Kn^ - 4.4 x 10^, Sc^ = 0.73) to large 

particles (g=10^®, R^ = 157 pm, Kn^ = 4.4 x 10"*^, Sc^ = 4.6 x 10^). In this 
connection it should be mentioned that particles of about 1 ym diam (see Fig. 
3 and Table I), which lie within the ” transition regime," correspond to the 
size ranges: 

10"° < g < 10"" 

3.4 X 10~^ cm < R, < 7.3 X 10*“^ cm 
1 

2.04 > Kn. > 0.95 
1 

1.1 X 10^ < Sc. < 5.35 X 10^ 

1 

(ii) As expected physically, the diffusivity of large particles is 
small, increasing significantly as particles approach molecular dimensions. 

On the other hand, the magnitude of the Soret factor shows an opposite 
trend in this case. The physical implication^^ of this is better understood 
from Eq. (11) which gives the net mass flux. Since is always positive a 

negative (or temperature gradient) tends to reduce the flux. In this 
case, the usually large when is small, and vice versa, leads to the 
conclusion that the Soret flux will dominate the Brownian flux only for 
larger particles. 

Figure 7 reveals the pressure and temperature dependences of the 
transport coefficients. Due to the complicated functional dependence of 
and on Kn^ and Sc^, simple representations may not be possible for all 
particle sizes. It is clear, however, that over the range of temperatures 
1700 to 3500 K) and pressures 0.5 to 1.0 atm) of interest in silicon 
reactors, retains approximately the same characteristics as described by 
the CE expression (i.e., T^^^/p), at least for particle diameters up to about 
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TABLE I 

COMPARISON OF VARIOUS SIZE RELATED VARIABLES FOR Si-Ar 




(P = 

1 atm, 

T = 2200 K) 



g 


R. 

1 

(cm) 

Kn. 

1 

Sc. 

1 

1.00 


1.58 

(-08) 

4.39 (03) 

7.2 

(-01) 

1.00 

(01) 

3.39 

(-08) 

2.04 (03) 

2.4 


1.00 

(02) 

7.31 

(-08) 

9.50 (02) 

7.8 


1.00 

(03) 

1.57 

(-07) 

4.41 (02) 

2.8 

(01) 

1.00 

(04) 

3.39 

(-07) 

2.04 (02) 

1.1 

02 

1.00 

(05) 

7.31 

(-07) 

9.50 (01) 

5.2 

(02) 

1.00 

(06) 

1.57 

(-06) 

4.41 (01) 

2.3 

(03) 

1.00 

(07) 

3.39 

(-06) 

2.04 (01) 

1.0 

(04) 

1.00 

(08) 

7.31 

(-06) 

9.50 

4.9 

(04) 

1.00 

(09) 

1.57 

(-05) 

4.41 

2.3 

(05) 

1.00 

(10) 

3.39 

(-05) 

2.04 

1.1 

(06) 

1.00 

(11) 

7.31 

(-05) 

9.50 (-01) 

5.3 

(06) 

1.00 

(12) 

1.57 

(-04) 

4.41 (-01) 

2.1 

(07) 

1.00 

(13) 

3.39 

(-04) 

2.04 (-01) 

6.8 

(07) 

1.00 

(14) 

7.31 

(-04) 

9.50 (-02) 

1.7 

(08) 

1.00 

(15) 

1.57 

(-03) 

4.41 (-02) 

4.2 

(08) 

1.00 

(16) 

3.39 

(-03) 

2.04 C-02) 

9.5 

(08) 

1.00 

(17) 

7.31 

(-03) 

9.50 (~03) 

2.1 

(09) 

1.00 

(18) 

1.57 

(-02) 

4.41 (-03) 

4.5 

(09) 
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TABLE II 


VALUES OF SOME INPUT AND COMPUTED CONSTANTS 


N3 


INPUT CONSTANTS (Si-Ar) 


p 

— X ci uITl 

P / T \ 

Gas Viscosity: — = { -zr— 

T 

= 2200 K 

Pr V / 

Om 

= 1.0 

Pr ~ 8.021 X 10~* poise 

“d 

= 1.0 

Tj. = 2000 K 

Ml 

= 28.09 1 , , 

\ gm/gm-mol 

0) = 0.35 


Me = 39.94 J 

otco = -0.10279 (CE theory estimate) 
N^ = 6.023 X 10^^ molecules /gm-mol 
k = 1.3805 ergs K"^ 




= 2.910 A 
(e^/k) = 3036.0 K 


COMPUTED CORRECTION FACTORS 


Oc = 3.542 A 
(Ec/k) = 93.3 K 
Ai/Ac = 350 
C 


^m 

Cd 

A 


= 1.0 
= 0.61557 
= 1.9182 


} 


(see Eq. (19)) 
(see Eq. (25)) 


Knc 


1.5 ) 

0.44 j 


(See Eq. (16) 


COMPUTED VALUES OF SOME VARIABLES 


Particle density: 

Pi = Pmp - Z(T - T^p)] 

V 

= 3.7498 

2 —1 
cm s 

^mp 

= 3 . 025 gm cm ^ 

p 

= 2.2116 

X 10*”^ gm cm^^ 

Z 

= 1.175 X lO"'* 

£ 

= 6.9572 

X 10“^ cm 

Imp 

= 1685 K 

Ac 

= 1.5472 

X 10~^ cals/ (cm-s-K) 
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1 ym. On the other hand, is nearly insensitive to temperature and pres- 
sure. Figure 7 is restricted to particle diameters up to about a micron 
since it is generally true that the Soret effect (thermophoresis) may not 
be as significant as inertial forces for very large particles. 

Figure 8 compares the case of a negative (e.g., Si vapor in Ar) 
versus one with a positive (e.g., Si vapor in H 2 ).* In order to increase 
the separation (or collection) efficiency of silicon, it is clear that 
should be positive so that the Soret flux is directed towards the ”cold” 
reactor walls. Since a multicomponent mixture of gases is actually present 
an ^’effective” Soret factor must be considered such that 



where x. is the mole fraction of species i and a., the Soret factor for species 
i diffusing in a carrier gas j . It is to be noted that the "effective*' Soret 
factor depends both on the relative concentrations of the various gaseous 
species present as well as on the magnitude and sign of their individual 
Soret factors. Typically, as in the Westinghouse process,^ for instance, 
one can expect a mixture of silicon in mainly sodium chloride vapor, argon, 
and hydrogen gases. The various mole fractions would be roughly 


'^Si 

= 0.08 


= 0.5 

^Ar 

= 0.12 

‘NaCl 

= 0.3 


If the effective carrier gas is a mixture of H 2 , Ar, and NaCl its average 
molecular weight would be 

M = h x.M. = 23.34 ^31) 

avg j J 3 


* 

Note, the value of for Si particles in both Ar and Ha is positive. 
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Since the molecular weight of Si vapor molecules is larger than this ^ 28) 

one would expect a positive value of to govern the Soret diffusion of 

silicon vapor in a mixture of the abovementioned gases. Indeed, using the 
above mole-fraction estimates together with the following approximate values 
of a (deduced from CE theory) applicable to Si vapor: 



0.7 

Si,NaCl 

- 0.3 

SijAr 

- 0.1 


one obtains s 0.27. Thus the net Soret flux of vapor should normally 

be directed towards the wall. Using this it was estimated that the Soret 
mass flux for Si vapor to the reaxtor wall may not exceed 15-20% of the cor- 
responding Fick flux. However, the sign of ex is sensitive to the pre- 

i,err 

vailing product concentrations. It is conceivable, therefore, that in certain 
portions, near the reactor walls, the concentration of H 2 may be locally small 
enough for to become negative and lead to a detrimental flux of Si 

vapor away from the walls. ^ To eliminate the occurrence of such a possibility 
we would recommend maintaining a Ha rich atmosphere near the reactor walls 
at all times - 

It is worth pointing out that the real problem, as mentioned 
earlier, involves (besides vapor) some silicon droplets that may be treated 
as "heavy molecules" (i.e., the submicron size class). In the context of the 
above discussion, pertaining to vapor, it is clear that the effective "molec- 
ular weight" of this "heavy-molecule-gas" will always be considerably greater 
than that of Si vapor (see Eq . (22))- Hence, one can, in fact, expect a huge 
(but favorable) Soret contribution to the overall mass flux for particles in 
the "transition regime." 

Finally, one needs to gain some idea of the relative importance of 
the three mechanisms of diffusion (i.e., Brox^ian, Soret, and turbulent) 
studied in this report. Bro^mian and turbulent diffusion are knomi to be 


Our more exact boundary layer calculations in the future will also be 
aimed at determining the extent of this possibility. 
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important for very small and very large particles, respectively. However, 
the importance of Soret diffusion for particles of about a micron diameter is 
easily seen from Fig. 7. Typically, ct^ is large ( =“ lO’) while is small 
( ** 10~*). Recalling that a^, the Soret factor, is related to a thermal dif- 
fusion coefficient D„ . (which is really the appropriate counterpart of D.), 

i , X X 

the following useful "thermal diffusion ratios" may be defined: 





a.Y. 

X X 


(1 - Y^) 


(32) 




T^ 


p(Di + Ep) 


T,i 


a.D.Y. (1 ~ Y.) 

1 1 X 1 


(33) 


The former compares the extent of Soret and Broraian diffusion, whereas the 
latter expresses the importance of Soret vs overall Fickian (i.e., Bro\mian + 
turbulent) diffusion. Since, as seen above, ^1, Y^ ^ 0.1, and 10 ^ 

(from Eq. (14)), close to the wall,* one finds: 


« 10 ^ 

« 10 

That is, Soret transport of silicon particles will predominate all the other 
forms of diffusion for particles about 1 ym diam, in a viscous sublayer region 
close to the reactor x^all. 

It is of interest to determine here the particle size ’h^indow” (see 
Fig, 3) for xdiich Soret diffusion dominates. An upper limit is approximately 
set by the physical requirement that particles exceeding about 1 ym diam xa/III 

19 . 

be controlled solely by their inertia. However, the condition k^ = 1 (see 
Eq. (32)) determines the lower size limit at which Soret diffusion becomes as 
important as Broxmian diffusion. It is interesting to note that this lower 
limit is sensitive to the prevailing Si concentration (or mass fraction) . 


Note that since Sp y^, the value of Sp drops rapidly as the wall is 
approached. For 0 < y“^ < 5, the variation in Cp corresponds to 0 < Sp < 
0.15 (see Table II). 
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Using — 0.1, one obtains ss 10, which (from Fig. 7) corresponds to a 
particle diameter of 0.01 ym, under typical reactor conditions. With = 10 
one can expect Soret diffusion to outweigh Bro\\mian diffusion, a condition 
which here would correspond to = 100 and a particle diam of 0.02 ym. Thus 
one can safely say that Soret diffusion (i.e., thermophoresis) will be the 
dominant transport mechanism, near the reactor walls, for silicon droplets in 
the size range: 0.01 ym to 1.0 ym diam. Note that this is a much wider size 
range than the one usually quoted in the literature on a somewhat ad hoc basis. 

In summary, the analysis just presented emphasizes the following 
aspects relevant to modeling the silicon separation processes within a reactor; 

(i) It is essential that the model describe the transport to the 
reactor walls of a whole size distribution of silicon particles, from vapor 
molecules to droplets of several microns in diameter. As a first step to 
determining the average mass flux due to a statistical particle size distribu- 
tion, the universally valid transport coefficient formulae, developed in this 
report, will be needed. 

(ii) Very small or large silicon droplets will be deposited on the 
reactor walls by Bro\smian and turbulent diffusion, respectively. However, 
the often ignored transport mechanism of Soret diffusion (or thermophoresis), 
due to large temperature gradients, will control the deposition of silicon 
di'oplets in the diameter range 0.01 ym to 1.0 ym, Xs^ithin a viscous sublayer 
region close to the x^?alls. 


32 


TP-392 


III, SILICON VAPOR DE POSITIOJ^ 

Using the GENMIX code we focus attention here on results of calculations 
which pertain to the mechanism of silicon vapor deposition and may thus help 
us improve the performance of a given silicon reactor, 

A, CROSS-STREAM PROFILES 

The structure of the developing, turbulent flow in the *'do^^mstream” 
section is characterized by the velocity, temperature, and silicon mass frac- 
tion profiles across the reactor tube radius. Figure 9 shows the variations 
in normalized (with respect to the centerline values) velocity, temperature, 
and concentration vs the normalized radial distance from, the pipe centerline. 
The operating conditions and pipe geometry were chosen to match the Westing- 
house reactor design.^ As expected from the prescribed mixing length dis- 
tribution discussed earlier, it is possible to discern at least three regions 
of distinctly different behavior. Near the pipe centerline (r/R - 0) the 
fully turbulent core flow tends to even out all gradients due to rapid mixing. 
On the other hand, near the wall (r/R = 1) the steepest gradients are estab- 
lished under the influence of the molecular dif fusivities for momentum, heat, 
and mass competing with residual (damped) curbulent diffusivities . Within 
this 'Viscous sublayer” the turbulent fluctuations decay rapidly (according 
to the Van Driest hypothesis) as the wall is approached. Thus, in the imme- 
diate vicinity of the wall, the final deposition of silicon may occur pri- 
marily by molecular mechanisms (i.e., a nearly laminar flow condition). 

Joining these two extremes is the outer turbulent (defect law) layer where 
transport occurs mainly due to diffusive turbulent fluctuations, characterized 
by a mixing length that varies linearly with distance from the wall. Note 
that in this region, unlike the viscous sublayer, there is negligible damping 
due to viscosity. Rather, the length and velocity scales of turbulence are 
smaller than in the core flow. 

Although in reality the abovementioned zones blend continuously into 
each other, it has been customary to view turbulent pipe flows in terms of at 
least three different layers (with respect to the velocity profile) : the 

laminar sublayer, the buffer layer, and the fully turbulent layer. Such a 
model is especially true of fully developed pipe flows. However, even Fig. 9 
(for a developing flow) reveals these zones. The profiles are nearly linear 
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next to the wall, implying that the laminar sublayer is a region of constant 
fluxes of momentum, heat, and mass. It is also interesting to note that the 
thicknesses of these three laminar sublayers are different, with the tempera- 
ture and mass fraction layer thickness being equal and larger than the velocity 
layer thickness. Such a behavior is dictated by the input values of the laminar 
Prandtl and Schmidt numbers. In this case Pr == Sc = 0,7 was used. On the 
other hand, the buffer layer region is also seen to be one of nearly constant 
fluxes, even though the fluxes here are less than those in the laminar sub- 
layf'r. Based on the prevailing wall shear stress, these two zones are expected 
to correspond to regions of the velocity profile that lie in the ranges 0 ^ y 
5 (laminar sublayer) and 5 y"^^ 30 (buffer layer). Both these zones fall 
within the viscous sublayer. Elsewhere in the boundary layer, the fluxes de- 
crease with distance away from the wall and the final boundary layer thick- 
nesses are nearly all equal; This is a consequence of the assumption that 
the effective Lewis and Prandtl numbers are nearly unity. 

B. STREAMWISE VAR IATIO N OF FLUXES 

Figure 10 shows the variation of the non-dimensional momentum, heat, 

and mass fluxes along the reactor length, X (non-dimensionalized with respect 

to the reactor diameter, D) . The flow and other conditions are identical to 

those used for Fig, 9. The following definitions of the Stanton numbers for 

heat (St.) and mass (St ) transfer, and the skin-friction’ coefficient (C ) , 
h m I 

were employed: 




w 


(34) 


- V 


St = 
m 


w 


b w 


(35) 




w 




(36) 
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X/D 

FIGURE 10 COMPUTED STREAMWISE VARIATION OF HEAT, MASS, AND 
MOMENTUM FLUXES TO THE RExiCTOR WALLS 

where the subscript "b" refers to bulk (average) quantities and "w" refers 

'\i 

to quantities evaluated at the reactor wall, h is the stagnation enthalpy 
which may be defined in terms of the temperature and flow velocity as 

h = T + Y~ ( 37 ) 

In the present analysis St^ = St^ due to the assumption of unity effective 
Lewis number. Furthermore, the closeness of the St and C^/2 curves in Fig. 10. 
suggests that under these conditions a Reynolds analogy of the type: 

St = C^/2 (38) 

might be valid within the reactor since the effective Prandtl (and Schmidt) 
numbers are close to unity and the pressure variation in the downstream direc- 
tion is small due to a low inlet velocity (in this case U , , = 20 m s~M. 

c,xnlet 

In the other cases studied, it was found that greater deviations from this 
Reynolds analogy occurred with increasing inlet velocity due to more signifi- 
cant pressure gradients. Table HI summarizes some of these results. It 
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TABLE III 

NON-DIMENSIONAL FLUXES 


U) 



U . - = 20 m 

c, ml et 

X/D 

s-‘ 

^f/2 

St, (= St ) 
h nr 

Nu„(= Na^)^ 

j:) 

O 

1 

1 

1 

1 

u lu^ 

- C D 

2. 5000E-03 

B. 3427E-02 

1. 1552E-01 

3_5594E+02 . . 

4. ,4Q1BE^D3 l . Q00BE±C0 

2. 5004E-01 

1 1321E-02 

1 2953E-02 

4. 0070E+01 

4. 4191E+03 

1. 0188E+00 

2. 7492E+00 

4. 714BE-03 

3. 1492E-03 

1. 5839E+01 

4. 3943E+03 

1. 0781E+C0 

5. 2479E+00 

4. 0373E-03 

4. 3114E-03 

1 31 19E+01 

4. 3468E+03 

1. 1201E+C0 

7. 7466E+00 

3. 645BE-03 

3. 8323E-03 

1. 1549E+01 

4. 304BE+03 

1. 1567E+00 ... 

1 0245E+01 

3. 4007E-03 

3. 3314E-03 

1. 0543E+01 

4. 2650E+03 

1. 1903E+00 

1. 2744E t-01 

3 2230E-03 

3. 3122E-03 

.J3. 801 lE+00 . 

_ 4_2273E+Q3 

JI_22J.9E+GD 

1. 5243E+01 

3. 0874E-03 

3. 1455E-03 

9. 22B0E+C0 

4. 1910E+-03 

1 . 2530E-I-C0 

1. 7742E+01 

2. 9807E-03 

2. 0147E-03 

a. 7699E+00 

4. 1558E-»03 

1.581CE+C0 

2 0241E+01 

2 B945E-03 

3 9094 E -03 

8 3941E+00 

4. 1217E+03 

1. 3091E+C0 

2. 2740E+01 

2. B237E-03 

2. B232E-03 

B. 0798E+C0 

4. 0BB5E+03 

l.a36SF-tC3 

2 5239E+01 

2 7646E-G3 

2. 7316E-03 

7. 8132E+00 

4. 0565E+03 

1 3633E+00 

2. 7738E+01 

2. 7146E-03 

2. 6916E-03 

_7_ 5B56E+00 . . 

_-4_Q261EeQ3 

l_3fi92F-k03 

3. 0237E+01 

2. 6726E-03 

2 6415E-03 

7 3936E+00 

3. 9985E+03 

1. 4131F+00 

3. 2736E+01 

2. 6369E-03 

2. 5995E-03 

7. 2345E+00 

3 975BE+03 

. 1.439SE+00 

3. 5235E+01 

2 6058E-03 

2 5638E-03 

7 1107E+00 

3. 9621E+03 

1. 4621E+00 

3. 7734E+01 

2. 5a07E-03 

2. 335ZE-03 

7. 0416E+C0 

3. 9679E1-03 

1 . 4a00£+C0 

4. 0233E+01 

2. 55‘?0E-03 

2 SlOOE-03 

7 0139E+00 

3 9920E+03 

1. 4926E+00 

4. 2732E+01 

2. 5397E-03 

2. 4875E-03 _ 

_7'-.D17.0E+0a _ 

__4_Q29BE+Q3 

__1_5QQ9E±C3l_ 

4. 5231E+01 

2. 5227E-03 

2 4A-.75F-03 

7. 0405E+00 

4. 0761E+03 

1. 5039E+00 

4. 7731E+01 

2. 507 5E-03 

2. 4496E-03 

7. 0768E+C0 

4. 1271E+03 

1. S037E+00 

5 0230E+01 

2. 4939E-03 

2. 4334E-03 

7 1205E+00 

4 lBOlE-t-03 

1 5101E+00 

5. 2729E+01 

2. 4B17E-03 

2 418BE-03 

7. 16B2E+00 

4. 2337E+C3 

1. S106L+C0 . 

5 522BE+01 

2. 4708E-03 

2 4033E-03 

7. 21B4E+C0 

4. 2871E+03 

1 3106E+C0 

5. 7727E+01 

2. 46092-03 

2. _930E-03 

._2_2703E+0Q 

_,4. 34Q2E-tQ3 

_l_3102E.-tGD 

6 OOOOE+01 

2 4r>22E-03 

2. 3e32E-03 

7 3196E+00 

4. 3876E+03 

1 5096E+CD 


Nu is the Nusselt number. 

Re^ is the bulk Reynolds number based on pips diameter. 

1) /U. is tlie ratio of centerline to bulk velocity, 
cl) • 
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TABLE III (Continued) 


U =50 

c ^ inlet 


X/D 

. 500QE-03 . 

. 5004E-01 
. 7491E+00 
. 2480E+00 
. 7471E+00 
1. 0246E+01 

1. 2745E+01.. ^ 

1. 5245E-^01 
1, 7744E+01 
2. 0243E+01 
2. 2743E+01 
2. 5242E+01 

2. 7741E-^01 
3. 0241E+01 
3. 2740E+01 
3. 5239E+01 
3, 7739E+01 
4. 0238E+01 
4. 2737E+01 

4. 5237E+01 

4. 7736E*^01 
5 0235E+01 

5. 2735E+01 
5 5234E+01 
5. 7733E+01 

6. OOOOE-^01 


s 


1 


C 


f/2 


St, (= St ) 
h nr 


Nu, (= Nu 
h m 




b 




3.Jb930E-02 
6. 7661E-03 
3. 2103E-03 
2. 73B1E-03 
2. 4685E-03 
2. 2B93E-03 
2^1596E-03, 
2. 0606E-03 
1. 9821E-03 
1. 91B2E-03 
1. B651E-03 
1. 8205E-03 
1.7a30E-03.„ 
1. 7510E-03 
1. 7234E-03 
1. 6995E-03 
1. 67B6E-03 
1. 6604E-03 
1. 6444E-03 
1. 6304E-03 
1. 6183E-03 
1. 6071E-03 
1. 5985E-03 
1. 591 lE-03 
1. 5B42E-03 . 
1. 5738E-03 


3. 2027E-Q2 
7. 9653E-03 
3. 6452E-03 
3. 046BE-03 
2. 7058E-03 
2. 4799E-03 
2. 3170E-03 
2. 1926E-03 
2. 0945E-03 
2 0146E-03 
1. 9493E-03 
1. 8931E-03 
1. B469E-:Q3 
1. 8077E-03 
1. 7741E-03 
1. 745PE-03 
1. 7202E-03 
1. 69a3£-03 
1. 6797E-03 
1. 6636E-03 
1. 6497E-03 
I 6376E-03 
1. 62B6E-03 
1 . 6206E-03 
1. 6129E-03 
1. 60/1E-03 


_4_QQ96E+02 

6. 1594E+01 
. 2. 7992E+01 
2. 3198E+01 
2. 0449E+01 
1. 8616E+01 
1.72 R3F+Q1 ■■ 

• 1. 625BE+01 
1. 543BE+01 
1. 4765E+01 
1. 4199E+01 
1. 3721E+01 

1- 3 313E-t-Ql 

1.2960E+01 • 

1. 2650E+01 
1. 2377E+01 
1. 2134E+01 
1. 1917E+01 

_1_1223E+Q1 

1. 1552E+G1 
1. 1413E+01 
1. 1354E+01 
1. 1356E+01 
1 1397E+01 
,-l..i44QE+01_. 
1. 1483E+01 


__L-lDiQE+04. 
1. 1047E->C4 
1. 0970E+04 
1. 0877E+04 
1. 0797E+04 
1. 0724E+04 
_1_0656E±Q4. 
1. 0592E+04 
L. 0530E+04 
1. 0470E+04 
1. 0411E>04 
1 , 0354E+04 
_1_0297£±Q4 
1. 0242E+04 
1. 01B6E+04 
1. 0131E+04 
1. 0077E+04 
1. 0023E+04 
.B-9hS3E.-f£>3 
9. 9200E+03 
9. 8831 E+03 
9. 904BE-I-03 
9. 9609E+03 
1. 0047E+04 
. i.0133Et.04 
1 . 0208E+04 


_l._0Qn3E+C0 
1. OlOSC+CO 
1-0493E+C0. 
I 0783E-e-00 
1. 1035E+00 
1. 1265E+00 
_l_JL479EtOQ 
1. 1683E+C0 
1. 1B77E+C3 
1. 2065E+00 
1. 2248E+00 
1. 2426E>C0 
_l_260ufcjH30 
1. 27/lE+GO 
.1. 2940E+G3 
1. 3107E+CO 
1. 3272E+C0 
1. 3436C+CO 
_L_3398E+nO 
1. 3758E-»-00 
1. 3909E+C0 
1. 4007E+C0 
1 . 4068E+00 
1. 4093E+C0 
_l_.41l3EtC0 
1. 4126E+00 


TABLE III (Continued) 


^c, inlet ‘ 

100 m 


X/l) 

^f/2 

St. (= St ) 
h m 

2. 5000E-03 

1. 9136E-02 

2. 7157E-02 

2. 5006E-01 

4. 3963E-03 

1. 7131E-03 

2. 7511E+00 

2. 0307E-03 

9. 9197E-04 

5. 2513E+00 

1. 8244E-03 

9. 0956E-04 

7. 7513E+00 

1. 6919E-03 

a. 3500E-04 

1. 0251E+01 

1. 598 5E- 3 

8 1594E-04 

1. 2751E+01 

1. 5281E-03 

7. 0576E-Q4 

1 5251E+01 

1 4724E-03 

7. 6174E-04 

. 7751E+01 

1. 4273E-03 

7. 4204E-04 

2. 0251E+01 

1. 3899E-03 

7. 2546E-04 

2. 2751E+01 

1. 3583E-03 

7. 1135E-04 

2. 5251E+01 

1. 3312E-03 

6. 9919E-04 

2. 7751E+01 

1. 3080E-03 

6. 8067E-O4 . 

3. 0251E+01 

1. 280OE-O3 

6. 7972E-04 

3. 2750E+01 

1. 2710E-03 

6. 7105E--O4 

3. 5250E+01 

1. 2561E-03 

6. 650/E-04 

3 7750E+01 

1 2432E-03 

6 3904E-04 

4 0250E+01 

1. 2319E-03 

6. 5370E-04 
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TABLE III (Continued) 
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TABLE III 


n =400 m s ^ 

. I nleL 


X/D 

^f/2 

St (= St ) 
n m 

2. 5000E-03 

1. 0918E-02 

__ 1.48358-02.. 

2. 4995E-01 

5. 1554E-03 
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5. 2473E+00 
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2 7605E-03 
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3 7727E+01 
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2 721 lE-03 
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should ba pointed out that *’bulk” values of properties differ from their 
corresponding centerline values (denoted by the subscript ”c ’) . This is 
clear from the variation of '^ith downstream distance. Typically, the 

centerline values monotonically exceed the bulk values in the region of de- 
veloping flow and become constant only when ^’fully developed* flow conditions 
are established.* 


C. SILICON COLLECTION 

The ratio of the mass of silicon deposited on the reactor walls, due 
to the boundary layer convective-diffusion processes, to the total mass of 
silicon entering a given reactor may be defined as the **collection efficiency.'* 
This ratio can readily be calculated by integrating the downstream variation 
of the silicon mass flux to the walls, over the reactor length of interest. 

Figure 11 shows the results of such a computation. It is seen that the silicon 
collection efficiency, for a given inlet flow velocity, increases as the reac- 
tor length increases. Initially, this increase is slower than in the aft re- 
gions (as expected from the mass flux variation shown in Fig. 10). Furthermore, 
Pig* 11 reveals the important effect of changing flow residence times on the 
reactor *s performance. As expected, faster through-flows reduce the collection 
efficiency drastically since the flow time is much less than the time for 
effective silicon vapor diffusion to the reactor walls. 


It has been sho\^m^ that in a pipe, as the developing shear layers merge, 
exceeds (U^) fully developed P^ior to settling dox^m to a state of equality. 
Such details of shear layer interaction, however, cannot be predicted by the 
present mixing length model. 
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FIGURE 11 COMPUTED VARIATION OF SILICON VAPOR COLLECTION 
EFFICIENCY OF A REACTOR WITH LENGTH AND INLET VELOCITY 


To summarize the results of this section: 

(i) Silicon vapor deposition processes were analysed with due 

regard to the structure of the developing turbulent flow that 
prevails over most of a reactor’s length. 

(ii) The modified GENMIX code was used to provide, a eompu tatlona I 

capability that should be useful to the designer in optimizing a 
given reactor, with a minimum of experimentation. 

(iii) A reliable basis for assessing the silicon vapor collection 
efficiency of a given reactor has been established. 

Finally, it must be mentioned that Soret transport of silicon vapor was 
intentionally excluded from this analysis, since it is expected to be small 
in the absence of condensation- The effect, however, will become important 
when silicon droplets are considered. 
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IV. A MODEL FOR SILICON PARTICLE COLLECTION 

la this section we consider the introduction of particles into the exist- 
ing GENMIX model^ for the vapor collection efficiency of silicon reactors of 
the Westinghouse type. ^ The last section concerned the partial collection onto 
the reactor walls of silicon vapor present in the reactor ^s hot core flow (see 
Fig. 1 for a schematic of the physical situation). I^Hiile the deposition of 
vapor was attributed to the processes of convection, Pick diffusion, and turbu- 
lent (or eddy) transport to the reactor walls, the deposition of particles 
(silicon droplets) is expected to be controlled by additional mechanisms. In 
particular, one expects Brownian diffusion to govern the deposition of small 
particles (10“^-10“^ ym diam) while Soret transport will become more important 
for larger particles (10“^-1.0 ym diam), especially within the ^viscous sub- 
layer” region next to the reactor walls, where large temperature gradients 
exist. Away from the walls the particles will be transported entirely by 
diffusive turbulent velocity fluctuations. This later effect will always tend 
to drive the particles close to the walls. As was shown abo\re (Section II), 
near the walls the former two molecular mechanisms (i.e.. Brownian and Soret 
effects) will predominate. 

In order to describe the overall silicon collection efficiency (generally 
due to vapor and droplets) of a reactor, it is essential to model the deposi- 
tion behavior of a distribution of particles varying in size from molecular 
dimensions to about a micron diameter. In the model considered here it is 
assumed that the particle size distribution function (PSDF) for the silicon 
droplets entering the downstream section of the reactor for separation/collec- 
tion is kno\>m. The question, as before, is, ”What fraction of the silicon vapor 
and droplets entering the reactor is collected on the walls?” The model 
described here attempts to answer this question, with due regard to the various 
depositioi?. mechanisms mentioned earlier and the two-dimensional developing 
boundary layer flow within the reactor. In fact, as a result of the slowly 
developing boundary layers within typical reactors, most of the deposition 
occurs prior to the establishment of ”fully-developed” conditions (i.e., with 
no streamwise variations). Hence, the available simpler solutions which are 
valid for fully-developed, or partially developed (Graetz type)^^ flow cannot 
be used since such analyses only provide an asymptotic limit for the silicon 
mass transfer rate to the walls. Over the reactor’s separation/collection 
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length, the velocity, temperature, and concentration boundary layers are devel- 
oping simultaneously « Hence, numerical solution of the coupled system of 
boundary layer conservation equations is unavoidable in determining the cumula- 
tive silicon collection efficiency, which depends upon the variation of the 
silicon mass flux (at the wall) along the length of the reactor. 

This section outlines the general approach adopted in the numerical model- 
ing of the silicon separation/collection processes. Specific complications 
are isolated in order to show the nature of the problem and the methods used 
to overcome the difficulties encountered. Insufficient time was available to 
fully couple a code based on this model to GENMIX. It is hoped that in follow- 
on work this will be accomplished. 


A. PHYSICAL AND MATHEMATICAL ASPECTS 

The present two-phase flow model for describing the transport of 
silicon particles from the hot ( == 3500 K) core-flow of the reactor to the 
cooled walls ( ^ 1685 K) is based on the assumptions that: (i) the total parti- 

cle volume is considerably smaller than the volume of the surx‘Ounding gas phase 
so that particle-particle interaction is negligible; (ii) the concentration of 
the particles is small enough so that their presence does not alter the surround- 
ing gas flow field (i.e., velocity and temperature distributions); (iii) the 
particles are small enough that they may be treated as ’’heavy-molecules" in 
local equilibrium with the gas. 

These assumptions are usually justified in the downstream section of 
the silicon reactor being considered. The diluteness of the silicon droplets 
can be expected to increase with increasing distance from the droplet source, 
immediately following combustion. On the other hand, the nature of the pro- 
cesses leading to silicon droplet formation (gas-to-particle conversion via 
nucleation) ensures that the resulting particles will be of sub-micron sises. 


Under these assumptions, one can write the following particle mass 
conservation equation suitable to the developing, boundary layer flow within 
the reactor (axisymmetric) : 


pu 





pv 


3r 


I 

r 



( 39 ) 


where the radial mass flux is given by 
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Expressions for the Brownian dif f usivity , D^, and the dimensionless Soret fac- 
tor, x^^ere presented above. The turbulent diffusivity, is taken 

to be unaffected by the small particles being considered and hence given by the 

>v 

expressions suggested by Lin et al, ' Equation (39) is valid provided the 
par:- :.cles do not grow appreciably (via vapor collection on their surfaces) dur- 
ing passage through the particle diffusion layer. Such an assumption may be 
justifiable in view of the typical thinness of these layers, as seen later. 
Moreover, Eq. (39) applies to particles of a uniform size, present in the 
external stream, depositing by convection. Brownian diffusion, and thermo- 
phoresis on the reactor x^alls. In order to describe the deposition of particle 
size distributions, therefore, one needs to repeatedly solve Eq. (39) xo^ith the 
appropriate boundary condition on concentration at the diffusion layer *s edge, 
corresponding to each particle size class. That is, a PSDF expressing the 
variation of particle concentration vs size in the external stream is dis- 
crete zed into a finite number of independent size classes, x^ith particles in a 
given class being characterized by uniquely defined transport coefficients, 
and D^. The velocity and temperature fields within the turbulent boundary 
layers on the reactor walls have to be supplied to the particle diffusion equa- 
tion. In the present case, for this purpose, the entire velocity and tempera- 
ture profile history along the reactor tube length is independently computed 
using a modification of the GENMIX code^»^ for the finite-difference solution 
of coupled momentum and energy equations (including variable fluid properties 
and the damping of turbulence near the reactor x^/alls) . 

Even in the absence of turbulence, the effects of thermophoresis and 
Broxrmian diffusion can immensely complicate the analysis of silicon droplet 
deposition behavior x^ithin the reactor. Physically, one expects thermophoretic 
effects to be important in regions of the reactor flow field x<?here temperature 


This implies that the particles folloxv^ the turbulent eddy motions of the 
gas. Actually, particles near 1 ym diam might break away from the eddies 
due to their large inertia. Such inertial effects are excluded in the pres- 
ent diffusional model. 
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gradients are large (i.e,, near the walls). On the other hand, the larger 
particle size class, for which thermophoretic deposition is significant, corre- 
sponds to Brox<mian dif fusivities that are extremely small. There exists, there- 
fore, a thin sublayer region next to the walls where both Brownian diffusion 
and thermophoresis may be important for particle deposition. Outside this sub- 
layer, however, temperature gradients continue to remain strong so that thermo- 
phoresis remains important, and diffusion by BroTOian motion of particles 
becomes negligible. Thus one can consider the presence of an "inner layer," 
next to the walls, where thermophoretic transport of particles balances convec- 
tion and Bro\<mian diffusion. External to this is an "outer layer" wherein 
convective transport balances thermophoretic effects in conserving particle 
mass.* Although Eq. (39) is, in principle, general enough to describe both 
these layers simultaneously, from the standpoint of computational efficiency 
it is useful to consider the abovementioned division into different layers 
(due to their disparate thicknesses) . 

It has been established that, in analyzing the effects of thermo- 
phoresis, it is convenient to visualize the introduction into the boundary 
layer of additional pseudo "sink" and "suction" behavior. The former will 
tend to deplete particle concentrations above the reactor walls while the 
latter acts towards increasing them. The net result of these opposing effects 
of thermophoresis is to appreciably increase the collection rate for the 
larger particles. 

B . COMMENTS ON THE NUMERICAL SOLUTION _P_EOX.E J^^^^ 

A fex^ difficulties arise x^7ith the numerical solution of the particle 
transport problem, described in the previous section. They are briefly 
discussed in this section, together with strategies adopted to overcome them. 


* It is to be noted that these layers are not related to the "viscous sublayer" 
region next to the x<?alls where turbulent dissipation of eddies occurs. In 
fact, both the "inner" and "outer" thermophoretic layers may be embedded 
within the "viscous sublayer" region. Moreover, the division into "layers" 
adopted here is different from the asymptotic three-layer model implied by 
Refs. 28 and 29. 


47 


TP-392 


1. The Stiffness Problem 

l^en the characteristic time scales associated with convection, 
Broxmian diffusion, and thermophoresis in Eq. (39) become widely disparate, one 
is likely to encounter "stiffness** during the numerical solution process. The 
problem is manifested by unstable numerical solutions if the integration does 
proceed or by a failure of the integrator to march forward, along the 
reactor length. The former corresponds to a build-up of truncation error as 
the solution proceeds (often seen as a wildly fluctuating solution) while the 
latter shows up as an inability of the integrator to meet the specified local 
error tolerance. 

To overcome this problem, a "stiffness solver,'* based on Gear's 

3 1 

variable order, variable step-size, integration formulae is adopted. A 
method of lines approach is utilized, suitably discretizing the derivatives 
along the radial coordinate, in order to solve Eq. (39) as an initial-value, 
ordinary differential equation. The "solver" employs an implicit scheme for 
greater stability and to allow marching in larger step sizes. Both these 
features are crucial to coping with the stiffness problem. 

2. The Scaling Problem 

Due to the disparate thicknesses of the "inner" and "outer" thermo- 
phoretic layers, discussed earlier, one runs into resolution difficulties if 
the adopted radial grid spacing approaches the thickness of the inner layer, 
or exceeds it. 

In order to deal with this problem the following transformation of 
the radial coordinate is defined: 

r - r 

_ w (l,0\ 



where r, r^, and r^ are the radial locations corresponding to the point of 
interest, the reactor wall, and the edge of the particle diffusion layer. 

Note that the transformation restricts the domain of integration of Eq. (39) 
to the range 0:^f^l, regardless of particle size. It is, of course, neces- 
sary to specify the approximate thickness of the diffusion layer (i.e., 

5m = r,, - r ) . Since the thickness of the thermal boundary layer (i.e., 6^). 

tu w ni 

at a given x, is known via prior solution of the energy equation, one may 
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readily estimate 6^ using the knom Lewis number, for that particle size 
class, from the relation*: 


6 

m 




1/3 


(43) 


3 , Determination of the Variable Edge Boundary Condition 

At the above value of 6^ a concentration boundary condition needs to 
be specified. While the other boundary condition on Eq. (39) is simply speci- 
fied at the wall as - 0, the one at the "inner layer" edge must be obtained 
via a degenerate form of Eq. (39). As discussed earlier, the degenerate equa- 
tion used here ignores the diffusion term and further assumes that the trans- 
verse convection of particles in the outer layer is negligible compared to 
their axial convection. Thus one can write the "outer layer" equation as 


u 



K 


(44) 


where ^ ^ : 
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and analytically solve it to obtain the following expression for the concentra- 
tion at the outer edge of the "inner layer": 


Y. = 
1 



(46) 


where the subscript ’o* refers to some specified initial conditions (e.g., at 
the inlet to the silicon separation/collection section of the reactor) ,’t 


* A more cumbersome approach could have been to estimate 6^ by solving the 
boundary layer concentration integral equation. 

t A more exact alternative approach to modeling the outer layer would include 
the transverse convection term in Eq.(44) and, thus, require numerical solu- 
tion of differential equations for the concentration along particle 
traj ectories . 
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Equation (46) emphasizes that the particle trajectories immediately outside 
the "inner layer" may be nearly parallel to the reactor walls. This is to 
be expected when the opposing thermophoretic and hydrodynamic forces acting 
on the particle nearly balance each other. 

C. STATUS 

The mathematical development described in this section has been incorpo- 
rated into a code (MPDEU) which couples to the boundary layer (GENMIX) code. 
Unfortunately time has not allowed its use. 
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V, DESCRIPTION OF THE GENMIX^MPDEU CODE FOR SOLUTION 
OF THg SILICON TRANSPORT EQUATION 

A. MPDEU CODE STRUCTURE 

Figure 12 shows the structure of the newly developed code MPDEU for 
the description of silicon particle transport in the developing tubular flow 
encountered in the downstream reactor section. Since Eq. (39) is a nonlinear 
partial differential equation, the well-kno^vni PDECOL^^ is utilized advantageous- 
ly for its solution. In Appendix A the particle transport equation, Eq. (39), 
is recast into a form suitable for solution using PDECOL. The main advantage 
of using PDECOL, as opposed to a conventional finite difference procedure, is 
its superior computational efficiency. It also incorporates variable step size 
control and variable order integration procedures that help it cope with numer- 
ical "stiffness” problems. A listing of MPDEU is given in Appendix B. 
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< TCOEFF-^EDGEY 
PROFILE-*- EDDY 


UINIT 


FIGURE 12 REPRESENTATION OF THE STRUCTURE OF THE PARTICLE CODE, MPDEU 


The roles of the various subprograms that constitute MPDEU are 
described below: 

MAIN 

a. Reads list of input parameters. 

b. Reads from random file DATA the results generated by the boundary 
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layer code regarding the various dependent variable profiles, etc. at several 
x-stations along the reactor (ioC., streamwise history) stored in ARRAYo 

c. Generates input parameters for PDECOL and satellite subroutines. 

d. Checks on input through TOITE statements. 

e. Calls VALUES for giving results at specified stations. These results 
include the non-dimensional mass fraction profile and the dimensional profiles 
along the tube radius. 

PDECOL 

This is the partial differential equation solver based on the method of 
lines approach developed by Madsen and Sincovec.^^ It uses a finite element 
collocation scheme for discretization of the spatial cross stream variable as 
it integrates a set of ODE with x as independent variable. A full description 
of PDECOL is available in Ref. 32 o 

BNDRY 

This is a subroutine in the PDECOL package, needed for specifying the 
boundary conditions on the partial differential equation being solved (i.e., 
the particle equation, in this case). In this case one specifies the two 
boundary conditions — at the wall and at the species layer edge. The one at the 
wall specifies zero mass fraction there while the one at the species layer is 
specified according to the method given in Section IV. 

UINIT 

This is also a subroutine in the PDECOL package. It specifies the initial 
conditions on the partial differential equation being solved. These initial 
conditions must be at an initial r. point (which is not x = 0) , typically a 
small distance downstream from the leading edge of the boundary layer, in this 
case XD(1). 

F 

A subroutine in the PDECOL package which specifies the differential 
equation being solved. In this case one specifies the particle transport 
equation in the non-dimensional form: 
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3x 
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(47) 


where A, B, C, D, E are coefficients whose values are specified in Appendix A. 
The values of these coefficients are computed in subroutine TCOEFF while sub- 
routine PROFILE is called to establish the velocity and temperature profiles 
at the x-station of interest. These profiles are needed in the evaluation of 
the coefficients, at points along the tube radius. 

DERIVF 

A dummy subroutine that is required by PDECOL (see Ref. 32 for details). 
TCOEFF 

This subroutine computes the coefficients A, B, C, D, E of the particle 
equation at a given (x, y) location within the tube. Prior to evaluation of 
the coefficients, however, it is necessary to compute the particle Damkohler 
number (see Appendix A), the suction terms due to variable properties, and the 
Soret effect and the ^effective diffusivity” (sum of the Brownian and eddy 
diffusivities) . Using cross-stream profile information on the abovementioned 
quantities at a given x-station, one interpolates to establish values for a 
given Y station. This subroutine calls EDGEY to find the boundary layer edge 
concentration value. 
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EDGEY 

The first task this subroutine performs is to find the mass transfer 
boundary layer thickness (DELM) , using the specified Lewis number and thermal 
layer thickness. The latter half of this subroutine computes the particle mass 
fraction YM at the outer edge of the boundary layer. This subroutine codes the 
solution of the degenerate particle transport equation (Eq. (47)). 

EDDY 

This subroutine computes the eddy diffusivity based on the Lin et al^^ 
formulae which describe the variation of the diffusivity discretely in the 
laminar sublayer, the buffer layer, and the turbulent outer layer. 

PROFILE 

This subroutine computes the cross-stream profiles at a given x-station 
by linear interpolation from the two nearest neighbor x-station profiles. Then 
the temperature gradient profile is established by numerical differentiation 
of the temperature profile (using subroutine DGT3) . These quantities are 
utilised to establish the needed profiles of Damkohler number, suction velocity, 
effective diffusivity, and the pseudo-effective rate constant. 

VALUES 

A subroutine which is supplied with the PDECOL package for the computation 
of results, at specified gridpoints, using a spline interpolation technique. 

B. COUPLING OF MPDEU CODE WITH THE GENMIX X0D_£ 

The solutions to the gas-phase momentum and energy equations are 
found using the modified GEIMIX code (see Ref. 2). These solutions are stored 
on disk through the random I/O procedure available on the GDC computer. Writ- 
ing these solutions on a disk file is accomplished by the WRITMS command. Once 
stored in the ARRAY vector these solutions may be read (via the READMS command) 
into a subsequent, independent run of another program. In this case, MPDEU 
utilizes these solutions to compute the silicon particle mass concentration 
field within the reactor. Special coupling statements are, therefore, present 
in both programs. 
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c 



Two sets of input parameters are involved: those that are direct 

input to PDECOL, and (2) those that are needed for the additional computations 
performed by PDECOL (e.g., as in subroutines TCOEFF, PROFILE, EDGEY, AND EDDY). 
The first category parameters are fully discussed in Ref. 32 and need not be 
repeated. In the second category the following parameters need to be input: 


YMO ~ Edge mass fraction at initial x-station. 

YW - Wall mass fraction. 

RW ~ Radius of tube. 

RMO - Radial location of boundary layer edge at initial x~station. 

KPOL - Degree of interpolation pol3momial to be used in subprogram 
AITINT. 

NGRID - Number of radial gridpoints at which solution is stored by 
GENMIX. 

XLO - Reference length used to non-dimensionalize x coordinate. 

NXPT - Number of x-stations at which profiles are stored by GENMIX. 

PR - Prandtl number. 

SC - Schmidt number. 

AMUR - Reference viscosity. 

TEMPR “ Reference temperature corresponding to reference viscosity. 

OMEG - Temperature exponent in viscosity law: = (T/T^)'^. 

DIR - Reference diffusivity. 

AND - Temperature exponent in diffusivity law: D.D. = (T/T )^. 

X i^ r 

ALPHA - Soret factor for given particle size. 

RUNIV - Universal gas constant 

NGAS - Molecular weight of carrier gas. 

ITURB - Index to indicate if laminar or turbulent flow is involved. 

ITEST - Index to obtain extra output (i.e., additional quantities) in 
case one wishes to check/test the program. 

ICOUNT “ Index to count the number of times the species boundary layer 
thickness is computed in sub EDGEY. It is used to distinguish 
the computation of the edge concentration at the first x-station 
from subsequent computations. 

YI(I) - Input mass fraction profile at first x-station to serve as ini- 
tial condition in the solution of the particle equation. 
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In the version of GENMIX coinmonly used the input is provided through 
"arithmetic FORTRAN statements" embedded x;ritbin the body of the MAIN program. 
Due to the many different types of parameters involved (e.g., control indices, 
geometrical, physico-chemical, etc.) the authors of the program apparently 
found it more convenient to divide the MAIN program into several subsections 
(or chapters) , with the input specified appropriately within each subsection. 
While this approach is useful for self-instruction purposes, it appears that 
once GENMIX has been suitably modified for a given problem it would be desir- 
able to pull out these "aritlunetic statements" and replace them by READ state- 
ments. Then all the necessary input can be specified separately on data cards 
and one saves the additional costs of recompilation of the program every time 
one or more of the input parameters change their values (e.g., as in a paramet- 
ric study of the reactor performance) . 

This was done here too and the different READ statements provided 
appear under a series of input points (within the MAIN program) entitled 
INPUT 1, INPUT 2, etc. (A listing of the version of GENMIX used is given in 
Appendix C.) In what follows, these different input parameters are described: 

INPUT 1 

LESSON - Index used to aid self-instruction in the use of the program. 
Does not affect the results. 

KASE - Index for self-instruction. Does not affect the results. 

NSTAT - Number of x-steps between output of single variables (e.g., 
pressure, etc.). 

NPROF - Number of x-steps between output of array variables, or "pro- 
files" (e.g., temperature, etc.). 

NPLOT - Number of x-steps betx^een output of plot of velocity, tempera- 
ture profiles. 

ITEST - Index that controls "extra output" for program checking or 
testing purposes. 

lUTRAP - Controls action taken if negative velocities are calculated 
(see Ref. 2). 

MTEST - Number of x-stations from start of computation for which 
"extra output" (with ITEST 0) is desired. 
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NSMALL - Number of initial x-steps over which one wishes to safeguard 
against possible instability by retaining a small step size 
(1/10 of specified step size) ^ 

MASSAG - Index that controls the need to improve specified initial condi 
tions by making them continuous profiles. It was not found 
necessary in the present calculations to modify the uniform 
initial conditions specified at the tube entrance. 

IREC - Index that controls record number on which computations at a 

given x-station are to be stored (on disk file DATA) for later 
use by the MPDEU program. Set IREC = 0 initially. 

INPUT 2 

ICUNIF - Index that controls whether uniform or polynomial initial condi 
tions are going to be used (see Chapter 5 of program). 

NFRAC - A fraction less than unity used in case polynomial initial 

conditions are specified. It specifies the wall value of the 
dependent variable related to concentration (see Chapter 5 of 
program) . 

NWALL - Number of radial gridpoints from tube centerline up to which 
a uniform I.C. profile holds and beyond which a polynomial I.C. 
profile is specified. 

ICURR - Index used to alter the values of the dependent variables at 
gridpoints next to the boundaries. ICURR = 0 does not alter 
them from the original GENMIX formulation. 

POHLIC - Pohlhausen parameter used in the specifications of polynomial 

I.C.s in Chapter 6 of the program. 

ICBET - Index that controls whether it is necessary to improve I.C.s 
(see Chapter 6) at all. 

IMPROV - Index that is used in case improved (polynomial) I.C.s are 
needed (see Chapter 6). 

RERR - Error tolerance in the computation of the tube wall streamline 
radius (used in Chapter 6)* 

Note ; The abovementioned INPUT 2 parameters do not play an important role in 
the program if it is found that adequate solutions result from uniform initial 
conditions alone, as was the case in the reactor study described in Section III 
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INPUT 3 
N 

XULAST 

LASTEP 

XOUT 

XEND 

FRA 

ULIM 

TAN 

PEILIM 

AFAC 

AEXDLM 

KRAD 

CSALFA 

KIN 

KEX 

POWER 


- Number of cross-stream intervals (between radial gridpoints) 
to be used. 

- x-location at which integration is to terminate. 

- Number of steps in x-direction after which integration is to 
be stopped. 

- Distance x (meters) to end of outer confining duct. 

- Distance x (meters) to end of central pipe. 

- Ratio of step length to total grid width, DX = FRA*Y(NP3). 

- Constant used in calculating entrainment. 

- Tangent of semi-angle of duct. 

- Maximum allowable increase in the quantity (i|^^ - step. 

- Factor multiplying area change term. 

- Maximum desired value of dimensionless area excess/2ir (i.e., 
dimensionless AEX) . 

- Control index for plane or axisymmetric geometry. 

- Duct geometry parameter (i.e., cos ot) which indicates if walls 
are inclined or not. 

- Index to specify nature of I boundary. 

- Index to specify nature of E boundary. 

- Exponent used in prescribing the grid spacing in terms of the 
normalized stream function w. 


& 


INPUT 4 


0M(1) - 0 ) value at tube centerline. 

0M(NP3) - 03 value at tube wall. 

ISTEP - Number of integration steps completed. 

lEND - Number of steps to end of central pipe, otherwise lEND = 10000 

(the latter applies here). 

lAX - Number of steps to position where I boundary meets symmetry 

axis, otherwise lAX = 10000. 

lOUT - Number of steps to end of outer confining duct, otherwise 
TOUT = 10000. 

XU . - Upstream distance, x. 

DX - Integration step size, DX = XD - XU. 
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IFIN - Index that controls when calculations should terminate 

(IFIN ^ 0 terminates) . 

EWALL - Constant in logaritlimic velocity profile formula. 

INPUT 5 

MASSTR - Index to indicate if considered flox^ has uniform or variable 
composition. 

NOVEL - Index to dictate if velocity computation is required. 

NEQ - Number of differential equations to be solved. 

INPUT 6 (SI units are to be used here) 

GASCON - Universal gas constant. 

CFU - Specific heat of fuel. 

COX - Specific heat of oxidizer, 

CPR - Specific heat of product 

CMIX - Specific heat of mixture. 

WFU - Molecular weight of fuel. 

WOX - Molecular weight of oxidizer. 

WR - Molecular weight of product. 

\^IIX - Molecular weight of mixture. 

VISFU - Viscosity of fuel. 

VISOX - Viscosity of oxidizer. 

VISPR - Viscosity of product. 

VISMIX - Viscosity of mixture. 

HFU - Enthalpy of formation of fuel. 

AENER - Activation energy of gas divided by the universal gas constant. 
PREEXP - Pre-exponential factor in chemical reaction rate law^ : 

w- T = F p^ m- m expC-E/R"*!) 
fuel ^ fu ox ^ 

OXDFU - Molar ratio of oxidizer to fuel in reactants. 

MODEL - Index to show x>?hether flow is laminar or turbulent. 

INERT - Index showing whether the considered mixture is chemically 
reacting or not. 

PRO - Value of laminar Prandtl number (assignable in terms of a 

profile x^rith distance from the wall). 

PREFO - Value of turbulent "effective” Prandtl number (assignable in 
terms of a profile \d.th distance from the x<rall) . 
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Note' ? Several of the parameters in this section pertain to a chemically 
reacting system as envisioned in Version 4A of GENMIX, However, in the present 
study, we consider a chemically inert flow within the ”doxmstream” reactor tube 
section. As such, many of these input parameters have to be assigned suitable 
“dummy” values in order to reduce the available generality to the specific prob- 
lem of interest here. A typical set of input parameters (given later) illus- 
trates this. 

INPUT 7 
H 

AK 
ALMG 
FR 
UFAC 

INPUT 8 

UA, UB,1 
UC, UD J 

TA, TB,\ 

TC, TD J 
TWALL 

RA, RB,*^ 

RC, RD J “ 

INPUT 9 

F2A, F2B,'l 
F2C, F2D J 

OXA, OXB 
OXC, OXD J 

Note ; In INPUTS 8 and 9 reference is made to four incoming streams of 
reactants (A, B, C, and D) according to the general configuration treated in 
Ref. 5. However, for the present problem of a single stream, containing silicon 
in a "carrier gas" having H 2 , Ar, NaCl as the main constituents, one needs to again 
resort to certain "dummy" values for these input parameters in order to 


Recovery factor. 

Mixing length constant K (von Karman’s), see Ref. 5. 
Mixing length constant, X. 

Constant used to calculate minimum velocity gradient. 

Constant used to calculate lower limit to . 

'dy ' 


Velocities of A, B, C, D streams. 

Temperatures of A, B, C, D streams. 

Wall temperature. 

Inner radii of A, B, C, D streams viewed as tubes. 


- "Fuel" concentration in A, B, C, D streams. 

- "Oxidizer" concentration in A, B, C, D streams. 
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simulate the 
parameters) 

INPUT 10 

PRESS 

DPDDX 

TREF 


VISREF 

VISEXP 

VPB 

VPA 


INPUT 11 
NINT 

KPOL 

NPTS 




problem of interest here (see the typical values assigned to these 
shown later. 


- Pressureo 

- Pressure gradient, dp/dx. 

- Reference temperature used in viscosity law. 



- Reference viscosity of gas corresponding to TREF. 

- Temperature exponent in viscosity lax^7. 

- Pre-exponential constant appearing in Si vapor pressure law. 

- Constant appearing in the argument of the exponential in the 
Si vapor pressure law: 



Mgi(VPB) 


M . 

mix 


P 




- Number of intervals used for Simpson’s integration to 
establish collection efficiency. 

- Degree of polynomial used for interpolation in the above 
integration. 

- Number of points (tNINlfl) considered in abovementioned 
Simpson’s integration. 
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E, TYPICAL INPUT PARAMETER VALUES FOR THE S.rLICON„.REACTQR 
SIMULATION WITH GENMIX 


INPUT 1 

INPUT 

3 (cont 

’d) 

LESSON 

= 0 

PEILIM 

= 


0.01 

KASE 

= 11 

AFAC 

= 


0.1 

NSTAT 

= 100 

AEXDLM 



0.01 

NPROF 

= 100 

KRAD 



1 

NPLOT 

= 100 

CSALFA 

= 


1.0 

ITEST 

= 0 

KIN 

= 


3 

lUTRAP 

= 2 

KEX 

= 


1 

MTEST 

= 2 

POWER 

= 


1 

NSMALL 

- 100 





MAS SAG 

= 0 

INPUT 

4 


IREC 

0 







OM(l) 

— 


0.0 



OM(NPS) 



1.0 

INPUT 2 







ISTEP 

= 


0 

ICUNIF 

= 1 

lEND 

- 


0 

VJFRAC 

= 0.5 

TAX 

= 


0 

NWALL 

= 10 

TOUT 


10000 

ICURR 

= 0 

XU 

= 


0 

POHLIC 

= 7.0 

DX 

= 


0 

ICBET 

= 0 

IFIN 



0 

IMPROV 

= 0 

EWALL 



9 

RERR 

= 0.01 






INPUT 3 

INPUT 

5 


N 

= 77 

MASSTR 



1 

XULAST 

= 9.0 

NOVEL 

= 


2 

LASTEP 

= 10000 

NEQ 



3 

XOUT 

= 9.0 





XEND 

= 0.0 





FRA 

= 0.05 





ULIM 

= 0.05 





TAN 

= 0.0 
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INPUT 6 INPUT 8 (cont'd) 


GASCON 

== 

8314 

TA 

= 

3500 

CPU 

= 

1375 

TB 

= 

3500 

COX 


1375 

TC 

= 

3500 

CPR 

= 

1375 

TD 


3500 

CMIX 

= 

1375 

TWALL 

= 

1700 

WFU 

= 

28.09 

RA 

= 

0 

wox 


25.58 

RB 


0 

WR 

= 

2 

RC 


0 

\^IMIX 


25.58 

RD 

= 

0.075 

VISFU 

= 

0.00005 




VIS OX 

= 

0.00005 


INPUT 9 

VISPR 

= 

0.00005 

F2A 

_ 

0.1 

VISMIX 

= 

0.00005 

F2B 

= 

0.1 

HFU 

=5 

0 

F2C 

=r 

0.1 

AENER 


0 

F2D 


0.1 

PREEXP 

= 

0 

OXA 

— 

0 

OXDFU 

= 

0 

OXB 


0 

MODEL 

= 

2 

OXC 

= 

0 

INERT 

= 

1 

OXD 


0 

PRO 

= 

0.7 




PREFO 


0.86 


INPUT 

10 




PRESS 

= 101325 


INPUT / 







DPDDX 

= 

0 

H 

= 

0.9 

TREF 



AK 


0.435 






VISREF 

= 

0. 00005 

ALMG 

— 

0.09 

VISEXP 

— 

0,65 

FR 

= 

0.033 

VPB 


7.3166E10 

UFAC 


0.01 

VPA 

= 

46710 


INPUT 8 


INPUT 

11 

UA 

=: 

20 

NINT 

= 

49 

UB 

= 

20 

KPOL 

= 

2 

UC 

= 

20 

NPTS 

= 

100 

UD 


20 
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VI , THE CHEMPART CODE; MODELS 

One of the major accomplishments of this research program has been the 
development of the CHEMPART code, a much modified and augmented version of 
the AeroChem Low Altitude Plume Program (the LAPP code) . ^ The CHEMPART code 
is an axisymmetric marching (parabolic) model which, in addition to contain- 
ing the non-equilibrium finite rate chemistry routines of the LAPP code, 
includes a number of new routines which treat two-phase flow phenomena, in- 
cluding: (i) the exchange of momentum and energy between a particulate phase 

and the gas, (ii) the formation of this particulate phase, i.e., nucleation, 
and (iii) the growth of the particulate phase via both heterogeneous condensa- 
tion/evaporation and via agglomeration. Another major improvement included 
in the CHEMPART code has been the addition of routines which allow the treat- 
ment of enclosed flows (in addition to the free jet expanding flows originally 
treated by LAPP). 

In this section we will describe the numerous models which have been 
incorporated into this code. The structure of the code, the preparation of 
input data, the operation of the code, and the output information obtained 
from it will be discussed in following sections. Finally, we will discuss 
the application of the code to Na/SiCl^* silicon reactors. Throughout the 
discussion of the CHEMPART code, frequent reference will be made to actual 
FORTRAN variables and subroutines contained within the code. This will be 
done to enable readers of this report to modify the numerous models contained 
within the code to better fit their purposes and to allow users of the code 
to correct errors which they may find have been incorporated within it. In 
any large code such as CHEMPART, many options are available to the user with 
regard to both the preparation of input data and the particular phenomena 
which may be treated. Not all these options have yet been tested and, there- 
fore, errors undoubtedly exist within the code which will only be corrected 
as it is used over a wide range of problems. Our aim here is, therefore, to 
present the reader sufficient information to allow him some facility in find- 
ing his way through the code, to augment or improve those models which he 
finds inadequate, and to correct any errors which may exist within the code. 

To help the reader who wishes to modify the code. Appendix D gives an explana- 
tion of the important FORTRAN variables used in the code. Appendix E ccntains 
a listing of the code. 
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A. GOVERNING EQUATIONS 


For gas-phase systems in which no particulate phase exists, the following 
equations describe the free shear layer turbulent or laminar mixing of the co- 
flowing axisymmetric streams undergoing chemical reactions. For turbulent 
flows all properties are interpreted to be the mean (time-averaged) values. 

The effective (molecular + eddy) viscosity, y, is described by one of the 
phenou^enological expressions given in Section VI. C (below). Nomenclature 
is tabulated in the front of this report. 

Global Continuity 


(pu) + i ^ (pvr) - 0 


Conservation of Species 
3F. 


pu 


!!i * „ !!i . 1 X !!i\ 

9x 9r r 9r Y Pr 9r J 


+ w,- 


Conservation of Momentum 


9u , 9u 
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- - t i * (« I7) 


Conservation of Energy 
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--f 

r or 


NS ( 
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i (Le - 1) , h, X (p 


.1 


State 


P 


pW 

RT 


(48) 


(49) 


(50) 


(51) 


( 52 ) 


The conservation equations are solved subject to • following 
initial and boundary conditions: 


X 

r 


0 : 


0: 


u = u(r), = F^(r), T = T(r) 

A = II = ^ 

9r 3r 9r ^ 
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For free jets: 


r co; u ->• u^, F^ ■> (\)e» ^ 


(53) 


For snclosed flows: 


3F, 

r = r„: u 0.1 cr. s“\ F^ (F^^) or 0*, T 

w 


(54) 


For a free jet expansion, pressure may be specified in the axial direction 
according to, 

p = Co + cix + C2X^ + CaX^ (55) 


It is convenient to transform the equations into a streamline coordinate 
system and utilize the stream function, as the radial coordinate. The 
transformation from cylindrical (x, r) coordinates to streamline (x, 'V) coordi- 
nates (which automatically satisfies global continuity, Eq. (48)) is defined by: 


u/ IIL 

^ 9r 
^ 9x 

From Eqs. (56a) and (56b) we obtain. 


pur 
— pvr 


(56a) 

(56b) 



pur 

/a \ 


'F 

\d'v) 


Introducing Eqs. (57a) and (57b) into Eqs. (49), (50), and (51), gives: 


(57a) 

(57b) 


Species 



1 

a 

(—] 

Upur^ 

aF. 1 

X 

¥ 

an' 

[vr) 

'F 

aH< 


(58) 


(F.) is specified by the equilibrium vapor pressure of the species at the 
^ w 

wall temperature, T^, if the ith species is condensible. Otherwise the 
radial derivative is set to zero at the wall. 
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Momentum 


Energy 


8u ^ 
3x 


= - ^ + 1 8 f ppur^ 3u 1 

pu dx w 3’i' [ H' 34* j 


( 59 ) 
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The governing set of parabolic partial differential equations, (Eqs. 
(58), (59), and (60)), are rewritten in finite difference form and then 


solved using standard forward marching techniques 


3 j 3 3 


The chemistry terms 


(see below), in the species continuity equations are evaluated via implicit 

differences^; the diffusion terms in the species continuity equations are 
evaluated via explicit differences.^^ 

B. FINITE RATE CHEMISTRY MODEL 

1. Kinetics Model 

Five possible reaction types are currently included in the program: 
Reaction Type 

(1) A + B C 

(2) A + B + M C 


(3) 

(4) 

(5) 


A 4 B 
A 4- B 

A 4 M 


C 

C 


4 D 
4 M 
4 D 


4 E 


C 4 D 4 M 


In Reactions (2) and (5) M is an arbitrary third body. In this program, all 
species are assumed to have equal third-body efficiences; thus, in evaluating 

The net rates of production for all reactions are written 
below, in the form w^^) Rp^j) - RM^^\ 


( 1 ) 


w^^ ^ = k F F 

w Kj p 


F F 


C D 


K 
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(2) 


^f 

D® F F 
A B 

V' "c 

w 


W 

K WRT 
P 

(3) 

w^J) 

II 

F F 

^ A B 

K 

P 

(4) 

w^J) 

= ^f 

A B 

"fp' "c 

K RT 
P 

(5) 


1 

F 

^ A 

kfP^ F^FdRT 

w 


w 

K W 


To reduce round-off and truncation errors, forward and reverse rates 
and are computed separately for each reaction. All contributions to 

the molar rate of production of a given species are then computed and added 
algebraically to form w^ (I®OT(J)) in subroutine CHEM. 


The forward rate coefficient, k^, is expressed in the form 


and K is determined from, 
P ’ 


= AT ^ exp(B/RT) 


InK = - AG/RT 

P 


(61. 


(62) 


The rate coefficients are currently divided into seven types: 
Rate Coefficient Type 


(1) 

^f 

= 

A 

(2) 

^f 

= 

AT“^ 

(3) 

^f 

= 

AT- = 

(4) 

kf 

= 

1 

>0 

(5) 

^f 

= 

A e;xp(B/RT) 

(6) 

"f 

= 

AT“^ exp(B/RT) 

(7) 

^f 


AT-"/" 

2. 

The 'rmo dynamic Data 




The thermodynamic properties (specific heat, Gibbs free energy, 
and enthalpy) for each species are input to the program as 
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- h298^\ 

Cp , ^ ^ j and (h^ — h 29 s^)in tabular form as a function of tempera- 

ture. Linear interpolation is used to define thermodynamic properties at the 
local temperature. Such data can, for many species, be found in the JANAF 
Thermochemical Tables. Where no such source exists, these quantities should 
be estimated in the Lest way available. (It should be kept in mind that any 
reasonable estimate is almost certainly better than exclusion of a species 
from the reaction scheme.) 

C , TURBULENT EDDY VISCOSITY MODELS ^' 

The following eddy viscosity models^ are incorporated into the 
program via subroutine TURB. 

Model 1 (Ferri)^ ^ 

Initial region, 

= p e = aO. 00137 x |p u —pul (63a) 

t ‘ o o e e' 

Developed region, 

= p e = aK bx /2 Ip^u^ - p^u^| (63b) 

where bi /2 is the value of r where pu = (p u + p u )/2 and K is the eddy 
' o e e e 

viscosity coefficient, usually taken to be 0.025. § 

Model 2 (Ting/Libby) 

2 2 

= pe = aK 7 x/2|u^ - u^l p(^) (64) 


where 


= 2 


/; 


(Pjp) r' 


d r' 


(65) 


For enclosed flows only Model 6 (Donaldson-Gray) or the laminar model may 
be chosen. 


t 

§ 


Defined as region upstream of axial position where 


(u — u )/(u. ~ u ) = 0.95. 
o e j e 


Most of the models contain a numerical coefficient K which must be deter- 
mined empirically. The value K = 0.025, taken from Schlicting,^ ® has been 
Incorporated directly into the program. This can be changed by the program 
input data via an appropriate value for the additional constant, a, Eqs. 

(63) - (69). 
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and ri /2 is the value of n where u = (u^ + u^)/2. 
Model 3 


Initial region, 

VI = p e = aO. 00137 x p |u — u I 
t '^o' o e' 

Developed region, 

= aK ri /2 — u | 

where ri /2 is the value of r where u = (u ■“ u )/2. 

‘ o e 


Model 4 

Initial region, 

y^ = p e 

Developed region, 

\ = p ^ 

Model 5 (Ting /Libby) 

Initial region, 

= P £ 

Developed region, 

= p 0 

Model 6 (Donaldson/Gray 
Initial region, 

y^ = p e 


a0,00137 X p lu — u I 
e ’ o e' 

aK rx /2 P e|u^ - u^l 



aK(r 1 / 2 — r . ) p u — u I /2 
' in o e ' 


(66a) 


(66b) 


(67a) 


(67b) 


( 68 ) 


(64) 


(69 a) 


In the CHEMPART program, the specification of Model 5 means that Eq. (68) 
will be used in the Initial region and Model 2 (Eq. (64)) will be used in 
the developed region. This Is important for re-starting a problem in the 
developed region for x^rhich Model 5 was selected to run from x = 0. 
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For Ml/ 2 ^ 

1.2 

K = 0.0468 


M 1/2 [ - 0,0460 M 1/2 

+ 0.0256 M 1 / 2 ''] 

Mi/2 > 

1.2 

II 

0 

0 

M 

03 



(70) 

where Mi / 2 is 

the 

value of the Mach 

number where u = (u + 

0 

u ) /2 (i. e. , the 
e 

half radius). 

The 

speed of sound 

at 

the half radius, ai /2 

is expressed by, 


, , = f— — 

Cp - (R/W 1 / 2 ) W 1/2 

where W 1/2 and T 1/2 are evaluated at the half radius. In Eq. (69a), r^^ is 
the inner mixing zone radius and is defined as the value of r where (u^j— Ug.) / 
(Uj - Ug) = 0.95. 

Developed region, 

Vij. = p e = OK ri /2 p|u^ - u^|/2 (69b) 

Laminar Flow 

The temperature dependent viscosity of each gaseous species is included 
as part of the input to the code and the viscosity obtained from them. The 
individual species viscosities are in the form 


1/2 


(71) 



The viscosity is then 


i = 1 


,NS 


NS 

u = ) F . W y , 

g i—t X X 

i=l 


(72) 


(73) 


D. GAS/PARTICLE DYNAMICS 

Small particles will (i) follow the gas streamline as they are 
carried along and (ii) adjust rapidly to the temperature of the surrounding 
gas (although phase changes such as fusion may result in sizeable transfers 
of energy to or from the gas) . Thus little exchange of momentum or energy 
will occur (with the notable exception of phase change just mentioned) . Large 
particles, on the other hand, will slip in the gas and may find themselves at 
a temperature significantly different from the surrounding gas. Thus, to 
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treat flows in which heavy particle loadings result in large gas/particle 
interaction, it is necessary to modify the gas momentum and energy equations 
to include this interaction and to introduce the momentum and energy equations 
of these larger particles. The following terms (74) and (75) must thus be 
added to the momentum and energy equations ^ respectively, (Eqs. (59) and (60)) 
to account for this momentum and energy exchange. For cases examined to date, 
particles with Kn 2 are found to have substantially the same temperature and 
velocity of the gas. Therefore, in the present version of CHEMPART, small 
particles are considered to be those for which the Knudsen number, Kn, is 
greater than 2, 


Momentum 


NPG F f u m. 

-9/(2y) > ^ 

p r! 

J P J 


U — Ur 


(74) 


Energy 


NPG y„m. 

Pj P^ Si 

9/C2p)> 

-4— J p r . 

i "^P J 


“ “ “p j + 


c_ S. 


Tp - { 75 ) 


The indicated sums are over the particle mass classes. 

The species equations for each particle mass class can be written 
in a form similar to Eqs. (49) and (58), i.e., in (x, r) coordinates: 


Particle Species 

9 F 


p u. 


Pi 


p. 

1 


9x 


4* p V. 


Pi 


a F 


9r 


= 1/r 


3r 



+ w, 


P. 

1 


(76) 


The equations describing the energy (including radiation to the wall) and 
momentum of the particles themselves are the following: 


The condensation (evaporation) of particulate material releases (absorbs) 
latent heat. The heat is assumed to be transferred immediately to the gas 
phase. The resulting change in gas temperature is calculated in subroutine 
AGGLOM. 
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Particle Momentum 


3 



3x 


9 

2 





( 77 ) 


Particle Energy 



3 y 



2 


r . 

J 


C 

Pr 



3bct /t ^ 

r . p ( p . 
3 P \ J 



(78) 


Equations (76), (77), and (78) are solved explicitly for ^’r) * % ’ ’ 

Fj_ y± 

respectively. In CHEMPART, terms and Eqs. (74) - (78) are evaluated and used 
only for particle mass classes for which Kn^ ^ 2. Particle classes with larger 
Knudsen numbers are treated as gases having the gas temperature and velocity 
and contributing to the gas density, enthalpy, heat capacity, etc. For small 
particles (Kn. ^ 2), radiation is assumed to cool the gas, rather than the 
individual particles and a further term 


NPG 


4ire'^ (T^ 


T‘‘) r! 
w 2 



is added to the energy equation (Eq. (60)). 

The terms f„ and are defined as follows^’: For Re„ < 0.5, f = 

1.0/SCF.. For 0.5 < Re„ < 10, 

f = [1 + 0.376 Re + 0.225 (Re )^ ln(Re )]/SCF 

Pn Pn P^ Pi 


For Re > 10 


= 0.033 Re /SCF. 


Pj 


(79) 


(80) 


(81) 


The factor g is given by 

g = (1 + 0.336 Re°'^^ Pr°’^^)/[Kn, x (1 + 0.336 Re°'== Pr°‘®® + Pr/Kn )] 

^3 ^ ^3 (82) 
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The factor SCF. in Eqs • (80), (81) is the Stokes-Cunningham correction 
factor^*^ 

SCF. = 1 + 1,26 Kn. + 0.4 Kn. exp(-l.l/Kn.) (83) 

J J J J 

In Eq. (82), Pr is the laminar Prandtl number. The source terms Wp . in Eq. 

(76) will be discussed below in the sections on nucleation and particle growth. 

The calculations described in this section are performed in subroutines 
GPINT and PARTC in CHEMPART. 

E. PARTICLE SOURCES 

1. Nucleation 

In systems in which very refractory species with large surface energies 
are found and in which large supersaturation ratios occur, such as those of 
interest here, it is attractive to employ a nucleation model^^ in which a limited 
series of simple gas-phase addition reactions are written which result in a 
’^critical nucleus” above which growth to parcicle species with bulk properties 
is rapid. Such a model recommends itself because (i) it fits naturally into 
the gas-phase chemical kinetics format of the CHEMPART code, (ii) it avoids the 
use of "liquid drop" models which utilize bulk surface tensions which, for 
small n-mers (n 10-10^ atoms), are difficult to justify, and (iii) it avoids 
the crucial assumption of classical nucleation theory that the nucleation step 
is rate-limiting — an assumption which will often be false if very high super- 
saturation ratios are encountered. Such a model has recently been proposed 
and used to treat iron particle nucleation in shock-tube studies in which 
Fe(C0)s was thermally decomposed . ^ ^ 

The model used to treat particle nucleation is one which treats the 
series of Si n-mer reactions 

Si + Si + M Si + M (84) 

k n n+k 

The model allows the forward and reverse rate coefficients for Reaction (84) 
to be estimated for all n. Using these rate coefficients it is possible to 
obtain a gas-kinetically determined critical size (as opposed to one based on 
classical liquid drop concepts) and a nucleation rate. 

Examination of processes like Reaction (84) leads one to write a 
realistic detailed mechanism for them as 
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*^n->n+k 
+ Si^ 7 • — •* 



(85a) 


. M . M (85b, 

where Reaction (85a) leads to an excited cluster which may either decompose 
into its original constituents or may be stabilized by collision with another 
molecule via Reaction (85b) . The overall forward rate coefficient for the 
overall Reaction (84) is then given by the following equation for ii^ terras 

of the rate coefficients for Reaction (85a and b) . 


n+k 


= ^n->n+k 


1 + 


[M] 


^n-hk->n 


( 86 ) 


The forward rate for Reaction (84) is then 



n+k" 


^n+k 


(87) 


Following the model development of Bauer and Frurip, 

values of K k follows: k 

s, n*fk n-^n+k n+k-^n n->n-fk 

gas kinetic (each collision yields Si’ ) so that 

n*rtc 


^^we estimate the 
is assumed to be 


K _ ,1 ~ 'ir(r )^ V - , 

n->n+k n+k rel,n,k 


( 88 ) 


where the radius of a spherical molecule containing (n + k) atoms. 

(The bulk liquid density is used to obtain this radius.) The ratio 

K , ,, /fc is the equilibrium constant for Reaction (85a); its value is 

n-^+k n-Hc^n ^ 

given by 


^ .1 .1 
n-^n+k n+k*m 


r 


47'' / ,2 

\ dr r 

n+k 



^n,k \ 
ICbT j 


(89) 
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where U . is an interaction potential between the n-mer and the k-mer and 
n,k 

k„ is Boltzmann's constant. The stabilization rate coefficient, k , is also 
B s 

assumed gas kinetic and is given by 


K 

S 




rel,n+k,M 


(90) 


where r^ is the radius of the bulk gas molecules and is the 

average relative velocity of cluster and bulk gas molecules. Using Eqs. (88), 

(89), and (90) the overall forward rate coefficient is obtained 

once the interaction potential U , is found. U , is estimated by Frurip 

n , n , 

and Bauer by a Lennard-Jones potential 


where e , 
n,k 



is a function of the heats of formation of the n-mers 




} 


/2 


(91) 


(92) 


and 


AH’’ 

n 


n 


|l( 1 - n-^/") 





(93) 


where L is the latent heat of vaporization. 

The overall reverse rate, for Reaction (84) may be determined 

from the forward rate once the Gibbs free energies of formation, AG °5 of the 
n-mers are specified. The enthalpies of formation necessary to derive AG° 
are estimated using Eq. (93). The entropies associated with the formation 
of n-mers from Si atoms are also obtained from Ref. Aland are given by 


AS'^ 

n 



(94) 


for n^lO. For n:^10, the configuration terms, on the RHS of Eq. (94) 

(the second two terms) are replaced by values obtained from the graph (Fig. 7) 
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of Ref, 41. Using these estimates, the free energies are given by 

AG° = AH® - TAS® + n AG? (95) 

n n n 


The reverse rate coefficients are given by 


n+k-^n 


k ,1 HT exp 
n-^n+k 


AG® + AG® -- AG®^^ 
n k n+k 

RT 


(96) 


The rate coefficients of Table IV were calculated using this model. 

It can be seen that the forward rate coefficients rise fairly rapidly with 
increasing n while the reverse coefficients decrease rapidly, reflecting 
the greater stability of the excited clusters with increasing n. 

For n:^10, the reverse rate coefficients range up to — 10 ^ ml 
molecule“^ s”^. For a simple two-body reaction between neutral, small 
molecules such rate coefficients are unreasonably large (by > 2 orders of 
magnitude). Their large size here reflects the fact that, at the high 
temperatures involved, the excited clusters of Reaction (85a) usually under- 
go decomposition as they are formed, before collision with another molecule 
can stabilize them. Larger (n > 10) clusters live longer due to their 
ability to channel excitation energy into the many available bonds and thus 
for these, realistic two-body reverse rate coefficients are formed for the 
overall Reaction (84) . 

In order to use this model to obtain rates of particle formation it 
is necessary to define a critical size cluster analogous to the critical 
size nucleus obtained from liquid drop nucleation theory. In the kinetic 
nucleation model used here this critical size is determined by establishing 
where the "bottleneck” in the series of reactions represp;:i.’'d by Reaction 
(84) occurs. This is done by finding that n-mer for whic> ' ■ rates of 
reactions creating it are mosL nearly equal to the rates de^ l '■ /ing it. 

The critical size determined in this fashion is compared to that obtained 
from liquid drop theory in Table V. To determine the number of Si atoms 
in the liquid drop critical nucleus, the Kelvin equation is used for computing 
the critical radius, r*, 


* 


r 


RTPplnS 


(97) 
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TABLE IV 

COMPUTED RATE COEFFICIENTS FOR Si POLYMERIZATION 
Si + Si^ + M Si^_^j_ + M (M = Ar) 

T = 2400 K, P = 1 atm 

Forward Rate Coefficient Reverse Rate Coefficient 

n (ml^ molecule"^ s~^) (ml molecule"^ s~^) 


1 

4.45 

(-30)^ 

4.59 

(-7) 

2 

1.01 

(-29) 

4.46 

(-8) 

3 

1.53 

(-29) 

4.33 

(-8) 

4 

1.97 

(-29) 

1.57 

(-8) 

5 

2.35 

(-29) 

1.68 

(-8) 

6 

2.68 

(-29) 

3.90 

(-9) 

7 

2.97 

(-29) 

9.57 

(-10) 

8 

3.24 

(-29) 

7.32 

(-10) 

9 

3.48 

(-29) 

6.90 

(-10) 

10 

3.71 

(-29) 

2.46 

(-10) 

15 

4.68 

(-29) 

1.78 

(-10) 

20 

5.32 

(-29) 

9.18 

(-11) 

25 

6.16 

(-29) 

5.72 

(-11) 

30 

6.78 

(-29) 

3.97 

(-11) 

35 

7.36 

(-29) 

2.96 

(-11) 

40 

7.80 

(-29) 

2.31 

(-11) 

45 

8.42 

(-29) 

1.87 

(-11) 

50 

8.90 

(-29) 

1.55 

(-11) 


The notation A(-B) implies A x 10 
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TABLE V 

COMPARISON OF CRITICAL CLUSTER SIZES 


FOR GAS 

KINETIC 

AND LIQUID DROP 

NUCLEATION 

MODELS 




* 

r 

(nm) 

T 

(K) 

P(Si) 

(atm) 

Supersaturation 

Ratio 

Liquid 

Drop 

Model 

Gas 

Kinetics 

Model 

1800 

0.001 

325 

0.44 

1.19 

1800 

0.01 

3247 

0.32 

0.93 

1800 

0.05 

16230 

0.26 

0.79 

2400 

0.1 

45.7 

0.16 

1.11 


whtre Y is the surface tension, Wg^ the molecular weight of silicon, Pp 

the droplet density, and S the supersaturation ratio. The number of Si atoms, 

n*, is obtained using . , p ,r*, and Avogadro's number N , from 

b X p A 

n* = 4/3 IT (r*) ^pN^/Wg^ (98) 


The surface tension of liquid Si is obtained from Ref. 42. 


Y 


720 


1.67 - 0.67 (T/1685) 


(dyne cm~^) 


The density of liquid Si used is given in Ref. 43. 


(99) 


p = 3.0247 - 0.355 x 10"^(T - 273.15) (100) 

P 

3 U 

The supersaturation, S, is computed using JANAF thermo chemical data for 
silicon equilibrium vapor pressures for Si(g) , Si 2 (g), sad SiaCg) above Si(£) , 

From Table V it can be seen that the kinetic model used selects a critical 
size larger than the liquid drop model predicts. For these large supersaturations 
the liquid drop model predicts nuclei radii which are on the order of one Si atom 
radius (0.17 nm) if one uses the surface tension of bulk Si(£) in Eq. (97) 
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(e.g., Y = 681 dynes cm ^ at 1800 K) . The use of such a value for y is clearly 
meaningless here. Thus the selection of larger critical sizes by the present 
model is physically more realistic for cases where very large supersaturations 
are encountered. 

The nucleation model just described does, however, have several 
deficiencies. The first of these is incorporated in Eq. (94) which describes 
the entropy of formation of the n-mers from Si atoms. The expression given in 
Eq. (94) does not converge toward the bulk liquid entropy as n increases; there- 
fore free energies of the n-mers do not, in general, converge to the free energy 
of formation of the bulk liquid as n becomes large. In order to rectify this 
situation the free energy of large n-mers is found by interpolating between the 
free energy obtained by the model just described and the sum of the free energy 
of the bulk liquid and the surface free energy of the n-mer, assuming the bulk 
surface energy value, i.e., we interpolate between the model just described 
and the classical surface tension model for the free energy of the n-mers. 

Thus, AG° is given by the following- equation: 
n 

AG°/n = F (AH° - TAS° + nAGS)/n + (AG° + 4irN, (1 - F ) (101) 

n nnn 5,An n 

where the fraction given by 

F = (500/n)/(500/n + n/500) (102) 

n 

From Eos. (101) and (102) it can be seen for n-mers much smaller than 500 the 
model just described will specify the free energy of formation. For n-mers 
with many more than 500 constituents the free energy of formation will be 
specified by the classical expression involving the surface energy of the 
bulk liquid. Equations (101) and (102) insure that, as n ^ the free energy 
of formation will indeed approach that of the bulk liquid. 

Actual calculations performed in subroutine NTHERM of the CHEMPART 
code are a simplified version of the model just described. A major simplifica- 
tion involves an abbreviation of the reaction mechanism described in Eq . (85). 
The code as it is currently written considers reactions only between identical 
n-raers, i.e., k = n in Eq. (85). Furthermore, three-body mechanisms described 
by Eq. (85) are simplified so that, only for the initial reaction between two 
monomers, does the third body explicitly play a role. Reactions between 
higher n-mers are considered to be two-body reactions. Further simplication 
is made by employing the approximate formulas for the forward reaction rate 
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by Frurip and Bauer, i.e., if n < 100, = 1 x 10”^® (n + k) ^ ® , If 

n ^ 100, = 1 X 10"^^ (n + k) ^ ‘ ® , The three-body rate constant for the 

reaction between monomers is chosen to be 1 x 10~^^, The nucleation model 
therefore considers a series of reactions of the type given in Eq. (85), leading 
to the formation of a maximum n-mer size chosen by the user. Above this maxi- 
mum size, the condensing material is considered to have formed a bulk particu- 
late size. The smallest particle mass class in this particulate phase has a 
mass which is somewhat less than the mass of the largest n-mer formed in the 
nucleation model. In each step of the program, as the largest n-mer is formi^d 
via the nucleation model, it is converted into the smallest and next smallest 
particle using the following three equations: 


.N 

w = Q F p N, 
p 1 nmax A 


= (1 - Q) F p N. 

P 2 nmax A 


where 


Q “ " W /N.)/(m'i “Till) 

^ nmax' A'^ ^ 


(In these equations the subscript nmax refers to the number of constituent 
atoms in the largest n-mer.) 

2. Agglomeration 


(103) 

(104) 

(105) 


The particle agglomeration model incorporated into CHEMPART (sub- 
routine AGGLOM) considers (i) Bro\mian agglomeration, (ii) curbulence-enhariced 
agglomeration, and (iii) agglomeration via inertial capture. To formulate 

the model, particles are divided into mass classes {m.}. - where 

1 1=1, r 

m (106) 


The factor 5 should be as small as possible for accuracy; because of calcula- 
tion time considerations it will probably be set at 10 which results in a 
rather minor overestimation^^ of the rate of growth of the particles. The 
rate of formation of particles of mass m^ via coagulation is given by 


.A 


w 


p± 


r i 


y K. . F F 

p,. 

j=l ^ 

+ y S. , K. . R. . F F 

^ ij ij P. 

j=l ^ 



• i-1 



+ y s. , . K. . . 
1-1,1 1-1,1 

i=i 

(1-R. .)p2 F F 

11 Pi P^ 



(107) 
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The fraction R.., required to insure mass conservation, is given by 
3-3 


R. . == 

13 


• . 1 •• (ni. + m. ) 

i-H ^ i 1 

m. , - — m. 

x+l X 


(108) 


is a factor which insures correct counting when particles of the same 
mass collide 


1/2 


S, . = 

iJ 


i = J 


± ¥ j 


(109) 


The rate coefficients for agglomeration, K.., are sums of coefficients repre- 

ij 

senting the various operative mechanisms, i.e,. Brownian diffusion, turbulence- 
enhanced agglomeration, and agglomeration due to particle-particle impaction. 


K.. - K®' + k’:”® + k“ 

IJ IJ 13 13 


( 110 ) 


Brownian Agglomeration ^ 


- «,T(r^+r )(D^ + D ) 


r . + r . 
1 - X 


13 


13 


4(D. + D.) ) 
1 3 1 

(r, + r^)Gy| 


( 111 ) 


where is the particle diffusion coefficient 


D. 

X 




67T y r . 
g 1 


1 + 1,26 — + 0.4 — exp 

r . r . ^ 

X 1 


(- A:) 


( 112 ) 


is the mean relative particle velocity 


13 


m. m, 

1 3 


8k^T 


1 /2 


(113) 


is a particle mean free path 


6.. = (6t + SI) 

13 13 


2 \ 1 / 2 


6 . = 

X 


21/2 

48 


7T G. 
X 

D. r. 

X X 




2 3/2 


(114) 

(115) 
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The following simplified expression is employed for certain cases: 
large particles with Kn^, Kn^ 0.01 


For 


K. . « 

13 


2 

3 y„ 


• V ft > t) 


(116) 


For Kn. or Kn. > 10 

1 j 


IJ 


■iT(r. + r ,) G, . 

1 2 iJ 


(117) 


Turbulence-Enhanced Agglomerat ion^ 


+ (Eo/v)^^“ (r^ + r^)" 3.j 


(118) 


where Eo, the turbulent energy dissipation rate per unit mass, is modeled by 


Eo - u, 


( 6)7 


(119) 


and is tha (eddy viscosity * molecular (laminar) viscosity) . The factor 
is the probability that a small particle approaching a larger particle 
will actually collide with and stick to the larger particle* This factor is 
a function of the Stokes number, given by the formula 

3 .. = 1 + 0 . 92/St. .+ 0.278/(St. .)^ + 0.034/(St. .)® (120) 

Kj ij ij ij 

obtained by fitting the (Re “ 0) curve given in Ref. 47. 

Particle-Particle Impaction 

For particles with differing average velocities, a contribution to the 
coagulation rate exists which is given by 


k“ . ,(r. +rj) = |>.p^-Up^|8y 


( 121 ) 


3. Heterogeneous Condensation/Evaporation 

Particle growth caused by the condensation of gas-phase species on 
the particles and particle shrinkage due to evaporation from particle surfaces 
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is treated using the following model. For large particles (Kn^ << 1) the 
rate of increase of particle mass is governed by Eq. (122)'*® 


( 122 ) 


where is the rate of increase of mass of a particle in the ith mass class 
and is the Sherwood number for such a particle. (^si^eq i 

Si mass fraction at the ith mass class particle surface which is taken to be 
the Si equilibrium pressure at the particle temperature (which is mass class 
dependent and different, generally, from the gas temperature). The gas density 
is p and is the diffusion constant for Si. The Sherwood number is a 

function of the particle slip (relative to the mean flow), Reynolds (Re), and 
Schmidt (Sc) numbers and is taken to be^® 


<®VKn<=l = 2(1 + 0.3 Sc/'*) 


(123) 


For very small particles (Knudsen number >: 1) diffusion will be very 
fast over lengths comparable to the particle diameter and the rate of mass 
increase will be limited by the actual incorporation of material at the surface, 
i.e., for an accommodation coefficient of unity. 


m . 


= Trr 


i(' 


Si 


- "si>eq 


.0 


8 RT 




l/2 


7T W 


Si 


(124) 


where W is the molecular weight of Si. 

o X 

The Sherwood number in this case comes directly from gas kinetics 
theory and is given by 




1 

2 


8 RT 




TT W_ . 
Si 


l/2 


r./D.. 
1 Sx 


(125) 


In order to compute the rate of mass increase throughout the broad particle 
size range desired, we will interpolate between the two regimes, using 




and 


m. = 2TTr. 

X 


. /y_. - (Y_.) Ap D_.Sh. 
X I Sx Sx eq, 1 1 Sx x 


(126) 


(127) 
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Condensarion and evaporation lead to a redistribution of particles 
among the various mass classes. This is accomplished by the following 
algorithm: 


For condensation (m.. 



0 ) 


w /p — 


m. 


pj.i j-i 


/(m^ - - 


"pj ■ "j’ 


(128) 


For evaporation (m. , 0) 

— — - - - - - J--- -- — 

.E 


w / P = F 


p^ - "j-i> 


- F ni. .^/(m - m.) 

Pj+1 3*^1 3-^1 3 


(129) 


In Eqs. (128) and (129) the factors involving m. , m , m. _ , m. , m . , 

3 J+1 3~1 3 J“l» and 

are those needed to conserve particle number and total mass. They are 

similar in nature to the factors R., used in Eq. (107) and are easily obtained 

^3 . 

by setting the sums of the particle mass after deposition (m^ + m^dt) equal 
to fractional parts of a particle in the jth mass class and C.he (j4*l)th mass 
class (for condensation) or (j-l)th mass class (for evaporation). 

As in the other "particle chemistry" models described above, this 
model rigorously conserves particle number- The code corrects the concentra- 
tion of the condensing or evaporating vapor so that mass is also conserved. 

The latent heat liberated or absorbed by the change phase is assumed to be 
transformed immediately to the gas at each step of the program and the gas 
temperature is adjusted (explicitly) accordingly in subroutine AGGLOM- 


F, ENCLOSED FLOWS 

A major problem encountered when considering enclosed flows is the 

calculation of the rate of transport of mass, momentum, and energy to the 

walls containing the flow. In order to do this a reasonable estimation of 

the turbulent transport (the eddy viscosity) must be obtained. In CHEMPART 

(subroutine TURB) the following simple mixing length governs the eddy 

viscosity used. The mixing length is related to the effective turbulent 

viscosity, M , in the wall region via Eq. (130) 
w 

y = y + p (130) 

w g 9y 
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where is the molecular viscosity and y is the distance to the wall. 


y 


w 


- r 


(131) 


The mixing length, Ji, is a function of the distance from the wall as given by 


the following two equations 

I = 


sa. 


K y n 
w 


A y,/K 
w Z w 


i = 


(132) 

(133) 

(134) 


for y 

In Eq. (132) the factor n is the van Driest damping factor given by*^ 

n = 1 - exp(-ypu*/A*^p ) 

S 

The constants and A^ have values typically taken to be 0.44 and 0.99, 
respectively. The friction velocity, u*, is (r/p)^^^, where x is the momentum 
flux to the wall (the wall stress). The boundary layer thickness is y^. The 
empirical factor A has a value which is a weak function of u* and dp/dx and 
is typically ^ 26.^ 

In the upstream regions of the reactor both boundary layer development 
and the spreading of an axially centered jet within the reactor must be modeled. 
In these upstream regions the following equation is used to describe the 
effective turbulent viscosity: 

M = [y.(r - r) + p (r - ri/ 2 )]/(r - ri/ 2 ) (135) 

J w w w 

In this equation is computed via the jet mixing model described in 

Section IV. C, while p is the turbulent viscosity given by Eq. (130). These 

w 

are averaged with weighting factors proportional to the distance from the wall 
and from the jet radius, ri/ 2 . l^en the jet radius exceeds half the distance 
between the initial jet radius and the wall, p is set equal to p^. 

Using the effective viscosity, P, obtained with this model, the code 
computes the flux of condensible species, particulate matter, momentum, and 
energy across the outermost stream tube to obtain rates of pressure drop, 
mass, and energy deposition on the reactor walls. For particulate transport 
to the reactor walls a simple thermophoretic diffusion velocity is also com- 
puted. The reader is referred to the discussion of the modified GENMIX code 
(Section IV) for a rigorous description of calculations of this type. 
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The change of pressure with axial distance, dP/dx, is computed using 
substantially the method of Patankar and Spalding^^ in which the change in 
pressure for a step, dx, is computed using the momentum flux to the wall, the 
gain or loss cf material via wall condensation/evaporation, and the temperature 
change which occurred during the previous step. The pressure change is that 
required to maintain the cross-sectional area of the flow equal to the reactor 
cross-sectional area. Since, in a marching code such as CHEMPART, the pres- 
sure change, dP, must be computed using infomation obtained in the previous 
step (i.e., downstream information required to precisely calculate dP is not 
available), the flow area will not be the same as the reactor area, but in- 
stead continually "seeks** the reactor area. 

G. GRIDPOINT DISTRIBUTION 

For free jet expansions the radial gridpoint in i|;-space is uniform, i.e., 
the stream tubes defined by the gridpoints all carry the same mass flux. As 
the jet entrains gas from its edge, additional gridpoints are added. However, 
the code, at present, is limited to a maximum of 30 radial gridpoints. Hence, 
it is necessary, when this limit is reached, to redistribute the gridpoints 
throughout the jet and reduce the number of gridpoints. At present, when the 
maximum number of 30 gridpoints is reached, this redistribution is done so that 
the code proceeds with only 20 gridpoints which are then augmented as more mass 
is entrained. 

For enclosed flows the amount of mass flowing through the system is con- 
stant (except for that condensing or evaporating at the walls of the system). 
Thus, throughout the calculation, a constant number of gridpoints is used 
and each stream tube defined by these gridpoints carries a constant amount 
of flow. For enclosed flow calculations, however, the spacing between the 
gridpoints is not uniform. Accurate calculation describing a number of 
phenomena occurring near the reactor walls requires that the spacing there 
be much finer than is, in general, needed within the core flow. Thus, the 
spacing (in 4^-space) of the first three gridpoints encountered as one proceeds 
inward from the wall is only 3.125% of that between gridpoints in the core. 

The spacing between successive gridpoints doubles as one proceeds inward 
until the 6th gridpoint is reached. The spacing is then even from the 6th 
point to the axis. In this way sufficient definition is obtained near the 
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walls to allow the calculation of reasonable values for mass^ momentum, and 
energy fluxes to the walls. 


H. STEP SIZE 

Except for the first few steps, the step size used by the program is 
totally determined by mixing parameters, in particular the value of the eddy 
viscosity and by the grid spacing, i.e., by (av* 0*- Specifically step size, 
dx, is governed by Eq. (136) 


dx = minimum [FDL ^ 


(136) 


where J/DL is an input parameter generally taken to be 1 and the subscript i 
denotes that values are examined at each gridpoint. For the initial fifteen 
steps taken in the calculation, a chemical kinetic criterion for step size 
is used in which the step size is taken to be a multiple of the molecular mean 
free path. Specifically, 


^ = X X ^(INTSTP - 1) 


U37) 


where A is the average of the molecular mean free paths at the centerline and 
at the wall or edge and INTSTP is the Ivjtegration step number. 
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VII. THE CHEMPART CODE: STRUCTURE 

A. GENERAL CHARACTERISTICS 

CHEMPART is written in standard FORTRAN IV. A listing of the code is 
included as Appendix E of this report. A simplified schematic of the pro- 
gram flows is shown in Fig. 13. The MAIN routine of CHEMPART governs the 
program flow. Within it numerous variables are initialized, routines for 
inputing and outputing data are called, and some of the major governing 
equations, i.e., the particle conservation equations, the momentum equation, 
and the energy equation are set up and solved. The initial routine called 
by the MAIN program is subroutine READ IT. Within this routine the input 
information is read and many additional variables are initialized. The next 
routines called are subroutines INITL, RADIAL or NONRAD, which set up initial 
values for concentrations, temperatures, and velocities in the coordinate 
system used within the program. After this input and initialization stage, 
the next subroutine to be called is subroutine NTHERM. This routine computes 
the thermodynamic data required by the nucleation and condensation models 
contained in the program. After a final initialization step performed by 
another call to subroutine INITL, subroutine WALCON is called (if the flow 
is an enclosed flow) and in this routine the vapor pressures of the con- 
densing species at walls is computed. The program then enters a large loop 
from which it does not exit until the problem is solved. Each passage 
through this loop marks the taking of a step downstream in the flow. With- 
in this large loop one first computes the gas-phase density, the molecular 
viscosity, and the specific heat of the gases at each radial gridpoint across 
the flow. (There are MPSI of these gridpoints.) After converting between 
the actual radial coordinates and the stream or ij; (PSI) coordinate system 
actually used to solve the problem, the subroutine TURB is called. In this 
subroutine one or another of the zero-equation turbulence models chosen by 
the user is called and the effective turbulent (or laminar) viscosity (the 
”eddy” viscosity) is evaluated at each radial gridpoint. The next step is 
the calculation of the step size to be taken (DX) . The MAIN program then 
solves the momentum and energy equations at each gridpoint. The equations 
are set up and solved implicitly for all gridpoints simultaneously. Sub- 
routine GETH is used during this step to obtain the static enthalpy at each 
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INPUT: SUB.READIT 


SUB.INITL 


.IRAO.^ 


SUB.NONRAD 


SUB. RADIAL SUB YFIT 


^nucl:: 


:WALL^ 


SUB.TURS 


FIND NEW OX 


SUB.NTKERM 


SUB.WALCON 


SOLVE GOVERNING EOS: 

IF .PART. « TRUE USE SUB.GPINT, 
PARTC., AGGLOM 


X: X-t> OX 


X>N*PRNT 


loUTPUT: SUB. OUTPUT 


X>XMAX 


STOP 


FIGURE 13 CHEMPART FLOW CHART 
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gridpoint and subroutine TRIDAG is used to solve the tri-diagonal matrix in 
order to obtain the new velocities and temperatures which result from taking 
the downstream step. The program calls subroutine TCHEM next. In this 
routine the species equations including chemical reaction source terms are 
set up and solved implicitly. MAIN then calls subroutine GPINT in which the 
contributions to the change in gas velocity and gas temperature due to the 
presence of particles moving through the gas at different velocities and 
temperatures are calculated. Subroutine AGGLOM is then called to compute 
the rates of particle agglomeration and particle growth via heterogeneous 
condensation. The particle conservation equations are then solved explicitly 
within MAIN and particle cloud densities are obtained for each particle mass 
class. Finally, subroutine PARTC is called and new values of particle tempera- 
tures and velocities are obtained using the momentum and energy exchange 
quantities found in subroutine GPINT. Subroutine NEGCHK is then called and 
the new downstream values for velocity, temperature, species concentration, 
and particle densities are scanned to determine whether a stable solution is 
being calculated (i.e., all these quantities are checked to see if any have 
become negative). If instability occurs, the step is repeated with a new 
step size (DX) which is one-tenth that taken initially. Subsequent steps 
then are gradually increased to that determined by Eq. (136). The above 
calculations are performed at each gridpoint at the end of which, if the flow 
is an enclosed one, subroutine FLUX is called from MAIN. Within subroutine 
FLUX the flux of energy, momentum, and condensible species to walls is com- 
puted and the pressure drop due to transport of momentum to the walls is 
computed. After return to MAIN, initialization required for taking the next 
step is performed, the program returns to the beginning of the large loop 
referred to above, and the next step is taken. This process is continued 
until the program has marched to the maximum downstream distance (XMAX) 
specified by the user. 

B. PROGRAM INPUT 

The preparation of input data for CHEMPART is described in Appendix F. 

A sample input card deck is given in Appendix G. The overall structure of 
the input data is governed to a large extent by the user’s choice of values 
for the logical variables IRAD, IRADP, PART, NUCL, and WALL (see explanation 
of Card 2 in Appendix F) . If the user chooses IRAD = .TRUE., then a set of 
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radial gridpoints Y(I), and values for gas species and particle concentrations, 
velocities, and tetpperatures at these gridpoints must be specified. If the 
user chooses IRAD = .FALSE., then the user must specify only axial and edge 
values for species and particle concentrations, velocities, and temperatures. 
For the latter case the radial dependence of the initial values of the quanti- 
ties will be step functions where the axial values are taken for gridpoints 
up to the edge of the jet (RJ) and edge values are used for all gridpoints 
with radius greater than RJ. The logical variable PART specifies whether or 
not a particulate phase is (or may be formed) in the flow. If PART = .TRUE., 
then initial values for particle cloud density, particle velocity, and parti- 
cle temperature must be specified. Whether they are specified only on axis 
and at the edge or at points, Y(I), across the flow is determined by the 
selection of IRADP to be either .TRUE, or .FALSE., respectively. The variable 
NUCL, if taken to be true, means that the code will include the nucleation 
model. If YIUCL = .TRUE., then cards specifying the number of n-mers to be 
included in the nucleation mechanism and various thermodynamic properties 
of the bulk liquid must be included (see Appendix G) . Modeling of particle 
growth via condensation and deposition of vapor on reactor walls is also con- 
trolled by the NUCL specification. If NUCL == .FALSE., these processes will 
not be considered. If WALL = .TRUE., then the flow calculation will be done 
for an enclosed flow and cards must be added to the input deck which give the 
reactor radius and the reactor wall temperature. 

C. PROGRAM OUTPUT 

The program output, a sample of which comprises Appendix H, begins with 
an identification of the program and an identification of the run (TITLE (I), 

I = 1,...36). Next follows information supplied by subroutine INCUT. This 
information includes the input pressure, the jet nozzle radius, the reactor 
radius (if an enclosed flow is being considered) , and the initial and final 
axial coordinates (XINIT, XMAX) . The print increment (PRNT) and the minimum 
step size (DXJIIN) are then printed. The viscosity model being used is then 
given. Next follows a table which gives the temperature, velocity, and 
species mole fraction along the jet centerline and at the edge of the flow 
(at R = infinity if a free jet is being modeled and at a point just beyond 
the jet boundary if an enclosed flow is being modeled). If particles are to 
be included in the calculation (PART = .TRUE.), a table is given of the number 
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of particle mass classes, the mass of the particles in each class, the radius 
of the particle, the particle cloud density along the centerline, the particle 
temperature, and the particle velocity. Next the reaction mechanism is tabulated 
each reaction is listed along with its rate constant. This list includes nuclea- 
tion reactions which are constructed by the nucleation model in subroutine 
NTHER^I. If particles are being treated in an enclosed flow (if PART = .TRUE, 
and if WALL = .TRUE.), the particle Schmidt numbers at the wall are then tabu- 
lated. The particle mass for each mass class is given along with the particle 
radius and the particle Schmidt number. 

As integration proceeds, whenever the x coordinate becomes greater than 
an integral multiple of the print increment (PRNT) , the following information 
is printed using subroutine OUTPUT: First, the value of the x coordinate is 

printed and the run labeling information is repeated. Next the non-dimensional 
value of X, the integration step size (DX) of the integration step immediately 
preceding the call to subroutine OUTPUT is printed, the pressure is printed, 
and the smallest and largest steps taken between print stations are printed. 

Next non-dimensional valued for the jet radius (QQIOO) , the inner mixing zone 
radius (QQ200) , and the Mach number at the jet radius (QQ300) are printed. 

Also, the mixing rate coefficient (QQ400) is printed if turbulence model 6 is 
used. Then the Prandtl number on the centerline is given. Next the non- 
dimensional coordinates of each gridpoint, the velocity, temperature, density, 
Mach number, static enthalpy, viscosity (eddy), and the value of the stream 
coordinate ^ are given at each gridpoint. This is followed by a tabulation 
of the mole fractions of each species at each gridpoint. If particles are 
present in the flow, the particle cloud density, particle temperature, and 
particle velocity at each gridpoint are tht=r> given. Next a summary of important 
physical quantities integrated across the flow is given. This includes the 
total mass flow (tt^^), the mass flow of gases, and the mass flow of particulate 
matter. The rate of change of the pressure and the total pressure change from 
the beginning of the calculation are printed. The total static enthalpy of 
the flow is then given followed by the total stagnation enthalpy of the flow. 
Changes in the mass flows are then listed: for total flow, gas flow, and 

particle flow. Next, changes in static and stagnation enthalpies and, if an 
enclosed flow is being modeled, the reactor radius and the reactor wall tempera- 
ture are given, followed by values for the mass average temperature and the 
mass average velocity across the flow. The energy and the mass being deposited 
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per unit length at the wall are then printed. The boundary layer thickness 
(non-dimensional) is then given. In the above, non-dimensional lengths are 
obtained by dividing by the initial jet radius if free jet expansion is being 
modeled or by the reactor wall radius if enclosed flows are being treated. 

The evaluation of integrated or averaged quantities across the flow is 
performed by making polynomial fits (subroutine POLFIT) to quantities which 
must be averaged or integrated and integrating the polynomials using function 
FINTGL. In order to obtain changes which occur and mass flow and enthalpy 
flow, it is necessary to subtract the values of the appropriate integral from 
the value of that integral obtained after the first integration step. It 
should be recognized that, early in the flow, the taking of differences be- 
tween large, nearly equal numbers obtained by such procedures is apt to be 
inaccurate. For example, changes in mass flow of 0.1% or less due to (i) 
entrainment of gas in a free expanding jet, or (ii) deposition of particles 
or vapor on reactor walls, will not necessarily be accurately calculated by 
the procedure used. 

A number of other output messages and output information will occasionally 
be given. These include messages concerning the presence of particles in the 
flow, a repetition of the particle Schmidt number information after the 100th 
integration step, and messages from subroutine NEGCHK if negative velocities, 
temperatures, species concentrations, or particle cloud densities are en- 
countered. Other messages which may be encountered in the output tell the 
user that the program has reached the end of the jet inviscid core, or that 
the number of gridpoints has been cut by subroutine TUBCUT. Lastly, there 
exist a number of STOP commands which may halt the program if certain errors 
occur. These STOP commands are numbered and the computer day file may 
generally be used to locate the source of the error with the help of the 
program listing. 
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VIII. SiCl4/Na REACTION AND Si(&) FUME FORMATION : 
SAMPLE CHEMPART CALCULATIONS 


In order to illustrate the use of the CHEMPART code we describe here a 
number of calcua.^ ons performed to describe SiCl^*/Na flow reactors of the 
type similar to that being tested by Westinghouse. ^ Specifically, a series 
of calculations in which the input enthalpy is varied over a substantial range 
will be described. For low input enthalpy the heat transfer to the walls is 
found to reduce the temperature in the core flow sufficiently rapidly that 
particulate formation occurs. Radiation from the Si(£) fume then reduces the 
temperature rapidly leading to even more fume formation and an accelerated 
temperature decrease. Calculations of this type test nearly all the models 
currently incorporated in the CHEMPART code; in fact, only the routines 
dealing with translation of the initial values from r-space to i(;-space remain 
untested at this time. In all the runs to be discussed in this section, step 
function initial conditions were employed. 

In Table VI the important input parameters of the series of five runs 
of interest here are described. The table shows input parameters for four 
runs in which the input enthalpy is varied by assuming various degrees of 
dissociation of a constant amount of input hydrogen. In the fifth run, more 
hydrogen was added and total dissociation of the input hydrogen is assumed. 

In the first four runs an Si input (in the form of SiCl^) of 28 kg h"^ is 
assumed; the fifth run assumes a flow of 45 kg h“^. The H/Si input ratio is 
4 in the first four runs and 12.5 in the last run. It is assumed that the 
arc heaters used by Westinghouse will supply between 0.3 and 1.9 MW in these 
runs. With these energy inputs, and starting with essentially room tempera- 
ture reactants, adiabatic flame temperatures varying between 2880 K (for 
run 1) and 2750 K (for run 5) are encountered. (Run 5 represents a case in 
which the input parameters match most closely those which Westinghouse plans 
to use. The first four runs underestimate the amount of pox^er to be added 
by arc heaters and the amount of hydrogen to be used.) Appendix G gives a 
listing of the input deck for run 1 of Table VI. The sample output shown 
in Appendix H is taken from the output of this computer run. 

The results of this set of computer runs can be summarized with the use 
of Figs. 14-17. Figure 14 shows temperature profiles along a streamline 
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TABLE VI 


SAMPLE TEST CASE - INPUT PARA^IETSE > 
Reactor Diameter = 15 cm 


^ c 
Input 


Run 

Si 

Input 
kg h"^ 

a 

Reactant Ratios 
H/Si Na/Si Ar/Si 

% of 

Input H 2 
Dissociated 

Initial 
velocity 
m s”^ 

Initial 

temperature, 

K 

Static 
enthalpy , 
kcal s~^ 

Heat 

added, 

MW 

1 

28 

4.0 

4.0 

0.5 

0 

23.1 

2820 

29.1 

0.31 

2 

28 

4.0 

4.0 

0.5 

6.2 

23.6 

2820 

32.6 

0.32 

3 

28 

4.0 

4.0 

0.5 

31 

25.5 

2860 

47.4 

0.38 

4 

28 

4.0 

4.0 

0.5 

62 

27.9 

2940 

66.8 

0.46 

5 

45 

12.5 

4.1 

1.1 

100 

125. 

3360 

380 

1.90 

^Reactants in 

runs 1-4 

are introduced as 

SiClz* (g) , H.2 , 

H, and Ar at 

3700 K and Na(g) 

at 1400 K. 

In run 5 


reactants enter as SiClz*(g) at 500 K, Ar and H at 4500 K, and Na(g) at 1400 K. Na(g) is introduced as a 
jet, r. = 3.75 cm; SiCl^,, H, H 25 and Ar are introduced in an outer flow between the jet and the wall. 

f 

Mass-averaged across the (initially) unmixed flow. 

j 

^This enthalpy includes heats of formation. ^ 

^The heat added is calculated by (i) assuming the H 2 and Ar must be heated to their input temperatures from 
298 K, (ii) SiCli*(£) at 298 K must be vaporized and heated to its input temperature, and (iii) Na(£) is 
supplied at 500 K and must be vaporized and heated to 1400 K. 
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FIGURE 14 TEMPERATURE PROFILES 

r/D ^ 0.27; wall temperature = 1800 ¥ 
Curve labels refer to run numbers. 
Conditions are given in Table VI. 
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FIGURE 15 Si(g) PROFILES 
r/D = 0.27. 

Curve labels refer to run numbers. 
Conditions are given in Table VI. 
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FIGURE 17 DROPLET SIZE DISTRIBUTIONS 

x/D = 3.4; r/D = 0.27. 

Curve labels refer to run numbers. 
Conditions are given in Table VI. 
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just outside the initial jet radius. It is along this streamline that reaction 
between the 3.75 cm sodium jet and the outer SiCl 4 /H 2 /H/Ar flow is greatest. 

As can be seen from Fig. 14, the computer runs to be discussed extend only 60 
cm or so down the reactor. The reason for not continuing these runs to great- 
er distances was simply one of cost. CHEMPART is an expensive code to run 
and the essential features of the flow were revealed in runs during the first 
50 cm or so. The essential features of the temperature profiles shown in 
Fig. 14 are as follows: For all but the last case, there is a rapid drop in 

temperature from the 3700 K input temperature for the outer SiCl^* /H 2 /H/Ar 
flow as H and H 2 rapidly react with SiClz, to produce, principally, SiCl 2 and 
HCl. As sodium crosses the streamline being tracked, the reaction to yield 
NaCl(g) and Si(g) causes the temperature to rise; however, for runs 1 and 2 
by x/D (where D is the reactor diameter) an Si(£) fume is formed which 
radiates strongly to the walls and ic> principally responsible for the de- 
crease in temperature noted in these two low enthalpy runs.* In run 3 con- 
densation does not occur until x/D 1.5. For the higher enthalpy run 4, 

particulate formation is not evident until x/D w 3 and even at x/D « 4, 

there is not a large amount of Si(£) fume. However, even for run 4 an 
accelerating temperature decrease may be noticed due to the small amount of 
fume present between 3 ^ x/D ^ 4. For the fifth run, no fume is formed 
and the temperature shows no sign of decay over the distance examined. In 

Fig. 15 the concentration of gaseous Si is plotted for runs 1-4 along the 

same streamline used in Fig. 14. In runs 1 and 2 the decrease of Si(g) as 
Si (5*) forms is very noticeable. For run 3 this decrease occurs later and 
is less drastic and for run 4 the decrease due to particulate Si formation 
is just noticeable by x/D = 4. [Si(g)] for run 5 is not shown; it is 
essentially constant and equal to 0.051. Figure 16 shows the Si(S^) cloud 
mass density along the same streamline. The rapid formation of particulate 
material in runs 1 and 2 at x/D < 1 is notable. For these runs, at x/D = 3 
about 3/4 of the available silicon has been incorporated in the particulate 


The model currently assumes that radiation occurs through an optically 
thin droplet cloud. For runs 1-3, the cloud is definitely thick optically, 
with an (absorption + scattering) coefficient » 1 cm”\ Thus the tempera- 
ture drop will not actually be as pronounced as indicated. 
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silicon. For the higher enthalpy runs, the formation of the particulate 
phase is still rather rapidly evolving even at x/D = 4. 

Figure 17 shows droplet SiC^*) size distributions at the streamline used 
for the previous figures and at x/D =3.5. As would be expected from the 

results shown in Fig. 16, the particle density distribution calculated in 
runs 1 and 2 shows considerably more particles present at this point in the 
flow than do the higher temperature runs 3 and 4. However, a more interest- 
ing feature is that the particle sizes present in the latter two runs tend 
to be larger than in the lower temperature runs. The particle size dis- 
tribution one obtains is expected to be dependent on the competition be- 
tween the nucleation process which generates new particles and the condensa- 
tion of Si(g) on particles already formed. Thus, in the higher temperature 
runs, fewer particles are formed and condensation on these particles tends 
to lead to larger particle sizes. In other words, the longer one is able to 
maintain the flow at relatively low supersaturations, the larger the resultant 
particles. The droplet sizes being calculated are relatively large. Par- 
ticle diameters in all four cases are considerably greater than 0.1 pm, and 
in run 4 the peak particle size approaches 1 ym. Indeed, the results shown 
in Fig. 17 hold promise of pointing to conditions at which particles sub- 
stantially larger than 1 micron might be produced by operating very close 
to the Si dew point. It should be recognized, of course, that because the 
size distribution is quite obviously a function of the competing rates of 
nucleation and particle growth via condensation (on the time scales involved 
in Figs. 14-17 particle agglomeration plays no role), it is important that 
the nucleation process be modeled reasonably accurately. The simple model 
used in CHEMPART is believed to yield considerably more reasonable results 
than would be obtained using a classical nucleation model with its assump- 
tions of rate limitation by the nucleation step and the use of surface tension 
concepts. The rate constants chosen for the nucleation reactions are de- 
liberately picked to be unreasonably large so that errors which may occur 
will tend to yield more> smaller particles than would be calculated if 
smaller rate constants were chosen. Nonetheless, the nucleation rates 
computed are slow enough that particles approaching 1 micron in size are 
predicted. Perhaps coincidentally, experimental evidence from experiments 
using SiClz*/Na flames at AeroChem, and even for seemingly unrelated silicon 
production processes such as the Union Carbide Free Space Reactor process 
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(an SiHi, decomposition process) , do seem to lead to particles with sizes in 
the 0. 1-1.0 pm size range. Certainly, more work to determine the sensitivity 
of the nucleation model to changes in the nucleation reaction rate constants 
is needed. Also, experimental work which would yield values for rates of 
nucleation in SiCl<,/Na flames is required before truly reliable results can 
be computed. 
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IX. CONCLUSIONS AND RECOMMENDATIONS 
A, CODE DEVELOPMENT 

In the research program described in this report, two extremely useful 
tools have been developed for studying and predicting the performance of high- 
temperature, two-phase flow reactors. The first, the combined GEN^lIX-MPDEU 
code, is suited for economically examining the development of the flow as it 
passes through the reactor and the transport of heat, condensible vapor, and 
particulate matter through the boundary layer to the reactor walls. The sec- 
ond, the CHEMPART code, treats the detailed chemical kinetics, particle forma- 
tion, and particle growth processes needed to adequately describe processes 
as complex as those occurring in, e.g., the SiCli*/Na reactor of Ref. 1. These 
codes complement one another in the following ways: (i) GEl'IMIX-MPDEU is not 

capable of performing detailed chemical kinetics including nucleation and 
particle growth calculations. However, by neglecting these time-consuming 
operations, it can, using a more detailed grid, perform its boundary layer 
calculations with more precision and economy than can CHEMPART. (ii) CHEMPART, 
on the other hand, because it treats many more potentially important processes 
than does GENMIX-MPDEU will, particularly in the core flow of the reactor, 
produce a much more complete picture of what can be expected to occur. The 
price paid for this is that a coarser grid must be used, with an accompanying 
loss of accuracy near the walls and, of course, increased run time and expense - 

Development of a computer code the size of CHEMPART is necessarily a 
multi-year project. It was noted earlier that many options for running the 
code and many types of systems can be treated by the code. To date, only a 
small number of these possibilities have been explored. It would be anticipated 
that, should the code find reasonably wide use, it will be substantially modi- 
fied in the future. At the present time, a number of possibilities for further 
improvement in the code exist. Also, the GENMIX-MPDEU code is as yet untested 
in its current form. In view of these points, the following recommendations 
are made: 

1. A sensitivity analysis of the nucleation model in CHEMPART should be 
performed. This should be done both by varying n-mer reaction rate constants 
and by specifying different numbers of n-mers (i.e., varying the size of the 
n-mer at which one treats the cluster as a particle of bulk material rather 
than a large molecule) . 
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2, There are a number of widely available integration routines specifically 
built to handle the type of stiff equations that must be solved in CHEMPART. 

Some of these stiff equation solvers should be tried in place of the mixed 
implicit-explicit integration schemes now being used in the code. 

3. The GENMIX-MPDEU code needs to be exercised so that accurate rates 
of particle transport to walls can be obtained for use in checking results 
obtained with CHEMPART and, in its own right, as a useful predictive tool for 
systems such as that of Ref. 1. 

B. CODE APPLICATION 

During this program the codes have been used to perform rather preliminary 
estimates of the behavior of flows in SiCl^/Na reactors of the general design 
of that described in Ref. 1. The results of the calculations may be summarized 
by noting: 

1. In flows with temperatures above the dew point of Si, velocities of 
less than 50 m s~^ are required to allow collection of more than 50% of the 
Si on the reactor (8 m long, 15 cm diam) walls. 

2. The calculations described in Section VIII show that large input 
enthalpies are required to suppress the formation of Si(£) droplets. The 
formation of a particulate phase at the high temperatures ( > 2500 K) pre- 
vailing in such flows is of particular consequence since the intense radia- 
tion from the droplets to the reactor walls represents a significant (actually 
dominant) mechanism of heat transfer. This is not unimportant with regard to 
mass transport to the walls, since, in order for thermophoresis to be effective 
in allowing the particulate phase to penetrate the laminar sublayer at the 
wall, a larger temperature gradient is desirable. Rapid heat loss due to 
radiation will reduce this gradient. 

3. In order to insure that sufficient enthalpy is introduced to prevent 
Si(£) droplet formation, large quantities of atomic hydrogen must be supplied 
via, in the reactor of Ref. 1, arc heaters. This results, for a given SiCl^ 
input, in an increased flow rate which in turn may lead to a loss in collec- 
tion efficiency for velocities ^ 50 m s The alternative is to decrease 
the SiCl 4 input, i.e., to decrease the production rate of silicon. One there- 
fore concludes that conditions must be sought which will allow the reactor to 
operate with the minimum input enthalpy per kg of SiCl^ input (chiefly in the 
form of H atoms) which permits negligible fume formation at flow velocities 

<50 m s*“^ . 
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A series of calculations describing the Westinghouse reactor is therefore 
needed in which the enthalpy is varied over a wide range, the H/Si ratio is 
varied over a reasonable range of about 4 to about 20, and the Si input rate 
is varied over a range of ** 25 to « 200 kg h“^. Calculations over several 
meters of the flow tube are needed to (i) check for the onset of Si(£) fume 
formation and (ii) determine wall collection efficiencies (the rate of Si 
deposition at the wall) . GENMIX-MPDEU could be used to assure that wall 
deposition rates (for heat and mass) are being performed accurately and to 
scan a wider range of conditions economically. Such a series of runs should 
provide an excellent prediction of conditions under which such reactors may 
be operated with greatest silicon production rates at optimum collection 
efficiencies. Other useful information, such as the rates of heat transfer 
to walls, pressure drop, and composition throughout the reactor would also 
be obtained. 


X. NEW TECHNOLOGY 


No reportable items of new technology have been identified. 
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APPENDIX A 


DFFTNITTON OF THE PARTICLE TRANSPORT EQUATION 

The particle transport equation, in primitive variables, has the form: 

3Y, 3Y. 9^Y. r,'" 

u + (v — V ) -5^ = D ^ i — (M) 

3x s 3r eff P 


3r" 


where, for the axisymmetric flow geometry. 


V = - D + 

s r eff 


^ i {<■ “eft} “1 (f S) 
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r.'" = p K Y. (1 - Y.) 


(pseudo suction velocity 

(pseudo volumetric reaction rate) (A3) 


with the pseudo specific rate constant given by 
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After the equations above are non-dimensionalized according to the following 
scheme: 
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(A5b) 

(A5c) 

(A5d) 

(A5e) 


one obtains the following definitions for the five coefficients that appear in 
the particle transport equation: 
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where 




(A6a) 
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Appendixes B through H are published 
in a separate supplementary volume, DOE/JPL- 
954862-79/8A. 


